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ABSTRACT 


Satellite x-ray observations indicate that most solar flares 
produce a hot, 2: lo^ K, dense plasma high up m the sun's corona. 

We believe that this gas originates from lower regions of the atmos- 
phere, specifically the chromosphere, where it is heated by a flare 
and evaporated into the corona. Optical observations have shown that 
the post-flare plasma does not cool uniformly rather, small cold 
condensations form at the top of magnetic field arches while the bulk 
of the plasma remains at high temperature. Loop Prominence Systems 
are an example of this phenomenon. Dynamical processes connected with 
the flare event have been proposed as a mechanism for these types of 
active prominences. 

We have investigated the cooling of post-flare plasmas, and 
attempted to explain the formation of loop prominences as due to a 

5 7 

thermal instability. At temperatures between 10 to 10 K, the 
solar plasma is known to be unstable to thermal perturbations because 
the radiation losses of this gas increase with decreasing temperature. 
A simple perturbation analysis is insufficient to account for the 
evolution of the solar plasma, because observed differences of temper- 
ature and density exceed an order of magnitude. Additionally, the 
plasma is always in a non-equilibrium state, so that, we must solve 
numerically the equations of motion and heat transfer. 

We have developed a one-dimensional model for active loop 
prominences. Magnetic fields present m post-flare regions are strong 
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enough to dominate the plasma, hence, only motion and heat fluxes 
parallel to the field need to be considered. The relevant size scales 
and time scales are such that single-fluid MH.D equations are valid. 

We have included in the model the effects of gravity, the geomeLry of 
the field and conduction losses to the chromosphere. A computer code 
for the solution of our set of equations has been constructed Basically, 
we treat the system as an initial value problem (with certain boundary 
conditions at the chromosphere-corona transition region), and use a 
two-step time differencing scheme. 

Our calculations indicate that the thermal instability mechanism 
is, by itself, sufficient to account for the behavior of active loop 
prominences. Under certain conditions, initial perturbations in the 
temperature and density profiles of small amplitude 5 $) and large 

size scale 2 X 10^ cm.) can grow into condensations with temperature 
and density differences of over an order of“magnitude and size scales 

g 

of less than 10 cm. In agreement with observations, the conditions 
that must be satisfied are such that Loop Prominence Systems are likely 
to occur only in large flares. The velocities, densities, and lifetimes 
that we obtain for the loop material are also in agreement with 
observations . 

From our results we conclude that the non-uniform cooling of the 
post-flare corona can be understood as a direct consequence of the 
temperature and density dependence of the radiative losses from a 
high-temperature solar plasma. 
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1 . THERMAL INSTABILITY 


1.1 Introduction 

Parker (1953) first introduced the concept of thermal instability 

4 

to explain the formation of cold condensations, T < 10 K, in the hot 
solar corona, T > 10^ K (Tandberg-Hanssen, 1967 )* He examined the 
stability of a hydrogen plasma to a thermal perturbation and deduced 
that if the total energy losses minus energy gains decrease with T, 
then the plasma would be unstable. This is because a region slightly 
cooler than its surroundings will lose energy faster than the surround- 
ings and, therefore, temperature differences will increase. 

Parker also noted that there is a temperature range over which 
the radiation losses of an optically thin hydrogen plasma decrease 

with temperature. This is the temperature range over which the state 

4 5 

of ionization of hydrogen changes rapidly, 10 ^ T ^ 10 . The emission 

is dominated by line and recombination radiation for these temperatures. 
Above this range, the emission is due mainly to free-free transitions; 
hence , 


X « 


„ 2 T 1/2 


( 1 . 1 ) 


where X is the radiative loss rate (Karzas and Latter, I 96 I). Below 

4 

10 K, the mean kinetic energy per particle is insufficient to excite 
neutral hydrogen above its ground state, and the emission decreases 
exponentially with decreasing T 

Defouw (1970b) has calculated the radiation from a model hydrogen 
plasma. He concludes that hydrogen is thermally unstable above 17,500 K. 
Thomas and Athay ( 1961 ) find that the solar chromosphere is unstable above 
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12,000 K. Cox and Tucker ( 1969 ) (Figure 1.1), Cox and Daltabuit (I 97 I), 
Tucker and Koren (1971) and Raymond et al. ( 1976 ) have calculated the 
total radiation losses (-including heavy elements), for the solar coronal 
plasma. From Figure 1.1 we note that there is a region, 10^ d T ^ 10^, 
over which S decreases with T, hence, thermal instability is possible. 

The original treatment of thermal instability by Parker (1953) 
was faulty because he did not take into account the dynamics of the 
plasma. Field ( 1965 ) m his excellent paper on the subject thoroughly 
analyzed the linear growth of thermal perturbations in a fully-ionized 
plasma. He not only included the dynamics, but also analyzed m detail 
effects due to conduction, magnetic fields, gravity, expansion and rota- 
tion. Goldsmith (1970) and Defouw ( 1970 a, b,c) extended Field's results 
to include ionization effects and Hunter (1970) investigated thermal 
instability in a plasma with an arbitrary velocity profile. 

We will deal only with thermal instability in post-flare coronal 
loops. The important effects in this case-are due to conduction and 
magnetic fields. Since the field dominates the plasma, we need only 
consider the growth of instabilities in one dimension, along the field 
(Section 3.4 ). Field has derived the growth rates for plasma in a 
uniform magnetic field, however, solar fields are far from uniform 
(Harvey et al., 1972, Rust and Bar, 1973)* This means that a constant 
flux tube in the solar corona will have a large change in cross- 
sectional area along its length. Hence, we derive below the growth 
equation (dispersion relation) for thermal perturbation m a loop of 
varying area, A(s). From the dispersion relation we can obtain criteria 
for a post-flare loop to be susceptible to thermal instability. 
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Figure 1.1 TEMPERATURE DEPENDENCE 'OF THE RAD I AT IV! 
LOSSES FROM A LOW DENSITY SOLAR PLASMA (Cox and 
Tucker, 1969)* 




1.2 Linear Theory 


The basic equations for fully- ionized plasma in a flux tube are 
discussed m Section 3.4 and Appendix A. They are 


d_p + 
at 


1 a_(Apv) 
A ds 


= 0 


( 1 . 2 ) 


<L_(pv) i a_(Apv 2 ) 

at a 9s 



9 (ppv 2 + , 1 SA fp pv3 + PV “ 

aF' 2 ^ 1 ) + a as' 2 Y-1 


f 3T 

a s 


+ X = 0 


(1.3) 

(1.4) 


and 


P 


R pT 
M 1 


(1-5) 


, where we have used Field's ( 1965 ) notation, except that we define X as 
the energy loss rate per cm rather than per gm, i.e. 


X(our) = p . X(Field )- . 


( 1 . 6 ) 


Note that m these equations we have neglected gravity. 

For the equilibrium configuration we assume that the plasma is static 
and homogeneous, and that the loo^ is infinite. Hence, the equilibrium 
state is given by. 


p = P 0 , T = T q , v - 0 , and X( p Q , T Q ) = 0 . (1.7) 

Following Field, we take the perturbation to be of the form. 


P(s,t) = P Q + ex P (Tit + iks) . 


( 1 . 8 ) 


For simplicity we assume that the area is given by, 
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( 1 . 9 ^ 


A(s) = Aq exp (k^s) , k^ = const. 


Substituting the form (1.8) into equations (1.2) - (l.5)j we obtain 


to first order* 



where 


k 


2 

1 


k ( k - *- k A ) 


( 1 . 11 ) 


and 


, -3_£(P,T) s i_£(p,T) 

jb SP * i T dT 


( 1 . 12 ) 


The dispersion relation is obtained from the requirement that the 
determinant of the coefficient matrix m equation (1.10) vanishes. 


As pointed out by Field, the dispersion relation has three solutions 
* 

one mode corresponding to isobaric condensations and the other two 
corresponding to the growth of sound waves. In this paper we will 
only be interested in the growth of condensations. We note from 
equation (1.3) that if the speed of sound is large so that the acceleration 
term is negligible, then we can replace the term (HPq) m the second 
row of equation (1.10) by zero. This corresponds to the case of 

5 
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isobaric condensations. Therefore, to isolate cases, we expand the 


determinant of the coefficient matrix about the second row, yielding 


h 2 p| 


_J3 

(y-i)t 


+ 


x T 


+ K k_ 


%pk. 


n_. k 2 P / jel 

'-1)T k l \ T 


(v-i)T 


+ 


X T + Kk L 


= 0 , (1.13) 


where we have dropped the subscript "0". The first term in (1.13) is 
the contribution due to the acceleration term in the force equation. 
If this vanishes then the dispersion relation reduces to 








(1.14) 


For instability, we need the real part of T| to be positive. From 
equations (1.11) and (l.l4), R g (t|) >0 implies 

h ~ T £ P 14 ' ^ • 0-15) 


Equation (1.15) is identical to the one derived by Field (I 965 ). The 
addition of a variable area has no effect on the growth of condensations 
other than to add an oscillatory factor to the time dependence. If 
conduction is negligible, then equation ( 1 . 15 ) reduces to 


dX _ £. ^ O 

ST T SP 


( 1 . 16 ) 


This relation was originally derived by Weymann (i 960 ). Note that 

Parker's simpler condition, < 0, implies that the radiative loss 

function for free-free emission, equation (l.l), is stable; whereas, 

under condition (1.16) it is unstable. However, we find that m the 

* 

corona, whenever the temperature is high enough so that bremsstrahlung 
dominates the radiation losses, the conductivity is usually so large 
that condition ( 1 . 15 ) implies stability. 
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We note from equation (1.15) that the effect of conduction is to 
increase the thermal stability of the plasma. This is to be expected, 
since conduction always acts to reduce temperature differences. For 
given values of n and T, equation (1.15) implies that there exists a 
maximum wavenumber, k c , such that perturbations of larger wavenumber 
are damped by conduction. Letting 


X = 2TT/k c (1.17) 

and using the Cox and Tucker ( 1969 ) form for the cooling function, 

X = n 2 A(T) , (1.18) 


-Of 

■where A 02 T , and O' is of order unity, we obtain for \ c . 


X « 2tt /-*f- . (1.19) 

V 3n A 

Hence, one of the criteria that a perturbation must satisfy m order 
for it to grow is that 


X > X 


( 1 . 20 ) 


Equation (1.19) provides a lower limit on the wavelength of an 
unstable perturbation. Field noted that the dispersion relation (1.13) 
also implies an upper limit on the wavelength of an unstable density 
perturbation. This is because the growth rate of a density perturbation 
cannot be faster than the speed of sound travel time across the pertur- 
bation; 


T = 1/kc . 
s 


( 1 . 21 ) 
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Therefore, f imte-speed-of sound effects will decrease the growth rate 
if the cooling time. 


T 

c 


a £ 

2 £ 


( 1 . 22 ) 


is less than t . This implies another criterion on the wavelength of 
an unstable perturbation, 

\ < ~ £C . (1.23) 

nh 

Note, however, that this restriction applies only to the growth of 
density perturbations. Temperature perturbations of arbitrary wave- 
length will grow on time scales of the cooling time as long as condition 
(1.20) is satisfied. In fact, we calculate in Section 5.4 that temperature 
variations m the post-flare corona grow at a much faster rate than density 
variations due to the finite speed of sound. 
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2. POST-FLARE REGIONS 


2.1 Solar Prominences 

Along with sunspots and flares, prominences are one of the most 
common manifestations of solar activity. In general, a prominence 
is simply a region of cold, dense plasma that is visible m Ha and is 
located at least partially in the corona. Since the corona usually 
has a temperature in excess of 10^ K, there is over two orders of 
magnitude temperature difference between a prominence and its 
surroundings . 

Prominences have a very wide range of morphology, however, it is 
possible to separate them into two categories : quiescent prominences 
which are large structures (size scales up to 10 ^ cm), that 
appear stable for months and occur in quiet regions; and active promi- 
nences which are short-lived events (lifetimes less than hours), and 
are believed to be flare related. (A thorough and detailed description 
of the observed properties of solar prominences is given by Tandberg- 
Hanssen (1974) in his book on the subject.) In this paper we will be 
interested mainly m active loop prominences. We include coronal ram 
and knots (Tandberg-Hanssen, 1974, P* 5) as part of this phenomenon. 

Dodson ( 196 I) first proposed that loop prominence systems are 
flare related events. A careful investigation of the association 
between loop prominence systems and flares was made by Bruzek (1964a). 

He concludes that they occur only after great flares and especially after 
flares that emit energetic protons (Bruzek, 1964b; Svestka, 1968 ). 

Teske (1971) found that loop prominence systems radiate m soft x-rays 
indicating that, in common with post-flare regions, they contain a 
high temperature plasma, T ^ 10^ K. A similar observation was noted 
by Waldmeier (1973)* 9 



We believe that loop prominences are a direct result of thermal 
mutability in hot, dense post-flare regions. They occur preferentially 
m these regions of the corona because the coronal density is greatly 
enhanced by a flare. The important features of flare regions and of 
loop prominences are briefly described m the following sections. A 
detailed description of the flare corona, including some of the recent 
Sky lab data is available m the book by Svestka ( 1976 ). 

2.2 Description of the Post-Flare Corona 

2.2.1 Magnetic Field 

The geometric structure of prominences and other features m the 
corona is almost wholly determined by the magnetic field. This is 
because the field is strong enough to prevent any mass or heat flow 
perpendicular to field lines (Section 3*4). Hence, temperature and 
density gradients perpendicular to the field tend to be much larger 
than those parallel. In addition, mass motion will be directed along 
field lines. 

The observed geometry of the field is a loop or a series of loops, 
each loop is believed to represent a magnetic flux tube. This is true 
if a flare region is observed either in soft x-rays (Kahler et al., 

1975, Vorpahl et al., 1975)» or m Ho? as a loop prominence system. 

The magnetic field strength in loop prominences appears to be large, 
as would be expected for flare regions. Zirin ( 1961 ) measured fields 
of up to 200 G, and Hyder (1964) obtained a value of 60 - 80 G for 
the field in a loop prominence. Similar values were observed by Roy 
(1972) and Rust and Bar (1973). Tandberg-Hanssen (1974, p. 43) estimates 
that the typical field in a loop prominence is from 20 - 100 G. 
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Several observers have compared the coronal field lines with those 
of a potential field computed from observations of the photospheric 
field. Rust and Roy (1971), Roy (1972) and Rust and Bar (1973). They 
find that the field in the post-flare corona can be well approximated 
as being potential, especially at large heights. This result is m 
agreement with the observation that the field m loop prominence systems 
is static. Bruzek (1964a) finds that individual loops do not alter 
their shape during the prominence event. 

In summary, we find that the magnetic field m the post-flare 
corona is large, 50 G, and appears to be current-free. 

2.2.2 Temperature 

Observations (Thomas and Teske, 1971) imply that all flares emit 

soft x-rays. The spectra of the soft x-ray bursts indicate that they 

7 

are due to bremsstrahlung from a thermal plasma at temperatures ^ 10 K. 
Typically these bursts have rise times of order minutes and decay times 
of order tens of minutes. The peak m temperature is usually observed 
to occur before the peak in flux (Horan, 1971; Kahler et al., 1970), 
which implies that the density of the soft x-ray emitting plasma is 
rising while the temperature drops. 

The dominant cooling mechanism for the hot coronal plasma is believed 
to be conduction to the chromosphere, Culhane et al. (1970), Moore and 
Datlowe (1975). A related process is the evaporation of chromospheric 
material. Although this has not been seen directly, it has been proposed 
by several observers as the mechanism by which the density m the flare 
corona is enhanced (Neupert, 1968 ; Hudson and Ohki, 1972; and Neupert 
et al., 1974). 
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The value of the temperature in flare regions is generally estimated 

either from the continuum or from lines in the soft x-ray spectrum. 

These measurements indicate that the maximum temperature is usually 

(2-4) X lo 7 K, but that the flare corona is highly inhomogeneous (Neupert 

et al., 1974, Dere et al., 1974). This is not surprising since, if a 

flare region consists of several loops, adjacent loops are thermally 

insulated from one another by the magnetic field. In fact, observations 

g 

indicate that loops of temperature >10 K can be seen in the immediate 
vicinity of Ha loops, T ~ 10^ K, (Fisher, 1971 » 1974, McCabe, 1973). 

If conduction is indeed the dominant cooling process, then we can see 
that large temperature differences would naturally arise. Even if all 
the loops m a flare region were heated to the same initial temperature, 

small loops would cool before large ones since the conductive cooling 

_2 

rate varies as L , where L is the loop length (Section 5.2.2). 

2.2.3 Density 

The density in soft x-ray emitting regions is not as accurately 
known as the temperature. From the observation of certain lines of 
highly ionized elements, e.g. Fe XXV, one can safely conclude that 
temperatures of at least 10 K are present; however, there does not 
seem to be any model independent method of obtaining values for the 
density. Part of the difficulty is that, like the temperature, there 
is no reason to expect a unique density for a post-flare region. 

From estimates of the emission measure, values for the density of 
10^ - 10 11 cm ^ are usually given, e.g. Cheng and Widing (1975). 

However, Craig and Brown (1976) point out that temperature mhomo- 
geneities may lead to large errors m the deduced magnitudes of both 
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density and temperature. In general, one may say that from soft x-ray 

10 -3 12 -3 
observations, n is at least 10 cm and may be as high as 10 cm 

m flare coronal regions (Svestka, 1976, p. 138). 

In loop prominences visible m Hey the density may be more easily 

obtained because the plasma is usually optically thick m the hydrogen 

lines. From the Stark broadening of the Balmer lines, Hirayama (1972) 

11 12 -3 

obtains densities ranging from 10 - 10 cm . Jefferies and Orrall 

(I 96 I) use the scattering of photospheric light to calculate the electron 

11 -3 

density. They obtain a value of 2 X 10 cm m a loop prominence. 

In summary, we conclude that the density m Ha loops is measured to 
11 -3 

be >10 cm , while soft x-ray observations indicate that the density 

10 -3 

m a general post-flare region is > 10 cm . Since loop prominence 
systems are only associated with large flares, it may well be that the 
higher value is more accurate for these particular flare regions. 

2.2.4 Morphology of Loop Prominences 

The evolution of a loop prominence system has been studied m detail 

by Bruzek (1964a). He finds that these systems first appear in Ha several 

' 9 

minutes after flare onset as a low mound, heights < lCr cm. The system 
then evolves as a series of Ha loops appearing at successively higher 
heights reaching up to 1.5 X 10 ^ cm. The apparent rate of expansion of 
the system is of order a few km/ sec, however, individual loops do not 
expand. The lifetime of a typical loop is approximately 30 minutes and 
the complete system may last well over 12 hours. Bruzek has observed 
certain systems lasting up to 3 successive days. 

Individual loops generally appear to begin as a bright small knot 
at the top of the loop, although knots can sometimes form at the sides 
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as well (Bruzek and Kuperus, 1972). The knot then expands into a loop 
with material streaming down both legs. This part of the evolution is 
very similar to "coronal ram" (Tandberg-Hanssen, 197^- > p. 29). The 
velocity of the material at the loop base is of order 10 km/sec, m 
agreement with free-fall velocities. The lifetime of a knot is of 
order 10 minutes; Bruzek and Kuperus (1972) quote 3-10 minutes, while 
Tandberg-Hanssen (1974, p. 5) states 15 minutes. 

2.3 Theories of Loop Prominence Systems 

2.3.1 Compression-Condensation Model 

Although Parker (1953) originally proposed thermal instability as 
a mechanism for the formation of quiescent prominences, he discounted it 
as a possible explanation for active prominences because their lifetime 
is much smaller than the growth time of thermal perturbations at coronal 
densities, ^ 10^ cm In addition, the size scale implied by equation 
(1.20) for coronal densities is much larger^than active prominence size 
scales. Field's (I 965 ) more exact analysis confirmed this result. 

Kleczek (1958) pointed out that if the coronal plasma were sufficiently 
compressed by magnetic fields, then the radiative cooling time may be 
decreased enough to agree with loop prominence time scales. Hence, a 

/ 9 "3 

model based on an initial strong compression (from n « 10 cm to 
n « 10 cm “’) has been investigated by several authors Kleczek (I 958 ), 
Lust and Zirin (i 960 ), Shklovskn (1965 )j Olson and Lykoudis ( 1967)5 
and De (1973). 

Several objections can be raised to this model. Jefferies and 
Orrall ( 1965 b) pointed out that the amount of matter m a loop prominence 
system, ~ 10 ^ gms, is too large to condense out of the corona, even in 
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a coronal condensation . Olson and Lykoudis ( 1967)5 themselves agree \ 

Q 

that the compression model cannot explain loops larger than ~ 3 X lCr cm, 
while loops larger than 10^ cm have often been observed (Bruzek, 1964a). 
Teske (1971) argues that the compression model is not consistent with 
his observations of soft x-rays from loop prominence systems. Finally, 
this model requires large changes in the magnetic field in order to 
compress the plama. However, observations indicate that the field is 
static during a loop prominence event, m fact, several authors find 
that it is potential (Section 2.2.1). Therefore, we conclude from 
these arguments that the compression-condensation model is not a likely 
explanation for the formation of loop prominences. 

2.3.2 Fast Particle Model 

Jefferies and Orrall ( 1965 a, 1965b) have proposed a model for loop 

prominence systems m which the prominence material is supplied by fast 

8 

protons, Vp ~ 1CT cm/sec. They argue that the large mass required for 
a loop system is stored as 10 keV protons on large loops extending into 

t 

the outer corona. These particles are supposed to migrate down to 
lower field lines where they thermalize near the top of a loop to produce 
the bright knots that are observed. 

The major objection to this proposal is that it is difficult to 
imagine how the protons can migrate across the field, and m such a 
way as to reproduce the evolution of a loop prominence system. Particles 
stored in the outer corona would have to migrate first down to the 
lowest loops and then to higher loops. Unless some mechanism whereby 
this may be done is proposed, the fast particle model cannot be truly 
evaluated. 
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2.3*3 Flare Evaporation Model 

Goldsmith (1971) and Sturrock ( 1973 ) have proposed that the mass 
for a loop prominence system is due to chromospheric evaporation by 
flare heating. They argue that since soft x-ray observations indicate 
that the density m post-flare regions is greatly enhanced over usual 
coronal values, a loop prominence system may well be simply the final 
stage of the evolution of the dense flare corona. The appearance of 
loops at successively greater heights may be explained as due to the 
decrease in the conductive cooling rate with height. 

The idea of chromospheric evaporation has another attractive feature 
m that one does not have to invoke magnetic field compression or other 
exotic mechanisms in order to obtain large enough densities for thermal 
instability to be effective. Goldsmith has calculated the response of 
a coronal loop to an isobaric perturbation in its initial temperature 
and density profile. He finds that for Tq — 5 X 10^ K, and n^ = 10^ cm 
a 5$ amplitude perturbation of wavelength, 3 X 10^ cm does grow to 
produce temperature differences of over an order of magnitude. However, 
these large temperature differences last for only < 100 sec, which is 
small compared to the lifetime of knots, ~ 10 minutes. In addition, 
his calculations do not include effects due to gravity, conduction to 

the chromosphere, or variation of cross-sectional area of the loop. 

7 

Therefore, the initial cooling of the 10 K plasma, and its effect on 
the formation of condensations cannot be investigated and the dynamics 

4 

of the cool, 10 K, plasma cannot be calculated. Hence, this model cannot 
be directly compared with the observations of loop prominences. 

We believe that the general features of the evaporation model 
proposed by these authors are correct. In this model loop prominences 
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are a natural part of the cooling process and are not due to some 
unobserved phenomena. Therefore, m order to critically test the 
model, we calculate the evolution of a post-flare loop including all 
the effects described above, and compare our results with observations. 
The model of a loop prominence that we use for our calculation is 
described in detail in the following chapter. 
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3. 


MODEL OF A LOOP PROMINENCE 


S KXfi® 


3.1 Introduction 

In this chapter we describe a one-dimensional (only variations 
parallel to B are considered), single-fluid model of a loop prominence 
which we have analyzed numerically. The physical assumptions that have 
been made and their justification are discussed below 


3*2 Single-Fluid Assumption 

The temperatures and densities used m the model are those observed 
for loop systems, 

3 X 1CT K ^ T 5 10 7 K , ( 3 . 1 ) 

11 -3 

and n ~ 10 cm . We assume that the plasma can be described by single- 
fluid equations. This approximation is valid only if the characteristic 
size scales and time scales of the fluid are larger than the mean free 
path and the relaxation time respectively. Spitzer (1962) has calculated 
the mean free path of electrons and protons in a fully ionized hydrogen 
plasma and the electron-proton relaxation time, 

\ £ m lO^ T 2 n" 1 cm. , (3.2) 

and 

T as 20 T^ 2 n 1 sec. (3.3) 

ep 

The length and time scales over which significant changes can occur in 

a loop are the thermal instability size scale, X c , equation (I.I 9 ) and 

the cooling time t , equation (1.22). For the values of T m (3-1) , 

we find that X £ /X and T /t < < 1. 

f c ep c 
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3-3 Radiative Losses 


Since hydrogen comprises ~ 90$ of the ions in the coronal plasma, 
the effects of the heavier trace elements are negligible except for 
their contribution to the radiation loss rate at temperatures above 

Ij. Gj 

~ 3 X 10 R. At high temperatures, (S \ty k), the radiation cooling 
is due almost entirely to lines m heavy ions such as 0^ + . These losses 
can be represented as a heat sink term m the energy equation, thus, we 
approximate the loop material by the fluid equations for a pure hydrogen 
plasma but with a cooling function, £, which includes the radiative 
losses of the heavy elements. 

The form of £ used is that derived by Cox and Tucker ( 1969 )* In 
their calculations they assumed that the plasma is of cosmic abundances, 
is optically thin to all radiation, and is m statistical equilibrium 
* with the states of ionization determined by a balance of radiative and 
dielectric recombination with collisional ionization. All the radiative 
and ionization processes considered vary as the electron density and 
ion density, hence the cooling function can be written as: 

£(n,T) = n £ nA(T) erg cm ^ s ^ , (3.4) 

where A(T) has been tabulated by Cox and Tucker, Figure 1.1. 

k 

The assumptions above are valid for loop plasma above ~ 3 X 10 K 
where hydrogen is almost fully ionized, but they break down at temper- 
atures below this. Once a significant amount of neutral hydrogen 
forms, the plasma will not be optically thin to hydrogen lines. The 
mean free path of a line photon is given by, 

X = 1/r^ ct , (3.5) 


19 



-14 2 

where a is the cross-section for absorption. For Ly 01, a ’ 10 cm , 

3 11 -3 

hence the mean free path is ~ 10 cm at densities of 10 cm . This 
means that if a slab of loop plasma only 1 km. thick were to cool to 

4 2 

temperatures ~ 10 K, the optical depth of the slab would be ^ 10 

Clearly, condensation in loops cannot be optically thin m at least 

Ly a, hence , in order to follow the cooling at low temperatures in 

detail one must include radiation transfer effects. 

Another assumption that breaks down at low temperatures is that 

the state of ionization is m equilibrium. Consider a plasma that is 

g 

cooling from a very high temperature, > 10 K, so that initially all 
the hydrogen is ionized. Assuming that the Cox and Tucker conditions 
are applicable, it will eventually cool to 1.6 X 10^ K with percentage 
of ionization, n /n « 50$. Now in order to justify using equilibrium 
equations to determine the ionization state, we must have the cooling 
time much greater than the recombination time, otherwise the temperature 
will change before equilibrium can be established. However, with the 
cooling rate and the recombination coefficient given by Cox and Tucker,, 

4 * 

we find that at.T 1.6 X 10 K, 



( 3 . 6 ) 


hence the assumption of statistical equilibrium is invalid at these 
temperatures . 

The arguments above imply that m order to calculate the behavior 

* 

4 

of the plasma below the 3 X 10 K level, we would need to follow it on 
very small size and time scales and include radiation transfer. There- 
fore, mainly for ease of calculation, we assume in our model that the 
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temperature does not fall below 3 X lCr K. This assumption is not 

unjustified since we expect that for the physical conditions present 

in loops, the important heat loss mechanisms, conduction and radiation, 

ii 

essentially vanish by ~ 10 K. Radiation losses become small since 
the temperature is < 1/10 the threshold energy required to excite 
hydrogen from its ground state, and the number of free electrons is 
reduced due to recombinations. Also, optical depth effects result 
in a decrease of emission. Conduction becomes negligible at low 
temperatures because the coefficient of conduction for a fully ionized 
plasma varies as a high power of T (Spitzer, 1962 ) 


K = 10"^ T 5 ^ 2 ergs cm -1 K 1 sec" 1 


(3.7) 


By fixing a minimum temperature we do ignore the profile of the 

4 4 

plasma between 10 - 3 X 10 K. However, we are primarily interested 

m determining whether the thermal instability mechanism can produce 


condensations in loops. As far as the hot loop plasma, > 10 K, is con- 

4 4 

cerned, the difference between 10 and 3 X 10 K is negligible. The 
main error is that we effectively overestimate the pressure of a 
condensation, hence we can say that if condensations do form under 
our assumptions, they will certainly form m a more exact calculation 
and will be cooler and smaller. 


3-4 MHD Equations 

4 

Since we restrict the loop temperature to T ^ 3 X 10 K, the state 
of ionization n £ /n > 99$* and we can take the plasma to be fully ionized. 
The single-fluid equations for a fully ionized hydrogen plasma in the 
presence of a gravitational field and a magnetic field are well known. 
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cf. Bragmskii (1965)* They are listed below 
the equation of continuity. 


P + v * (pv) = 0 ; (3-8) 

the equations of motion, 

p (|^+v-V)v+V p + V- n=pg + ijXB ; (3.9) 

the heat equation, 

— * 

p (|^ + v • V) e + P V . v + 7 . q + (n ' V) V = -£ + j [E + ~ X B] , (3.10) 


and the equation of state. 


p = ^ kT . (3-11) 

H 

The quantity e in (3.10) is the internal energy per gram of the plasma, 
and m our case. 


e _ 3kT # 

"h 

The stress tensor II and the heat flux vector q can be broken up 
into two parts, one parallel and one perpendicular to the magnetic 
field. 


(3-12) 


q=K VT+K 7 T , 
H it II x i 


(3.13a) 


and 


4r*- 44 


II = ri w + Ti, w , 
If II k 1 ’ 

Sv. 


(3.13b) 


a 

where the elements of W are terms like — — . For our model the 
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components of q and II perpendicular to B are negligible compared to 

those parallel. This is because the diffusion coefficients for a 

particular direction, 1C and ri depend on the mean free path in that 

II » 1 II » x 

direction. Parallel to the field the relevant quantity is the collision 
mean free path, and perpendicular to the field the gyroradius is important. 
Braginskn (I 965 ) obtains. 


K 

__1L 

K 


( "g T col>‘ 


1 

•n. 


where t . is the collision time and <0 is the gyrofrequency. 
col g 


(3.14) 

Taking 


10 G, and n 


10 ^ cm 3 , yields. 


(to T .) 2 = 10 6 B 2 T 3 n " 2 = lO" 1 ^ T 3 . (3.15) 

g CO i. 

For temperatures T ^ 10 K, the ratio of parallel to perpendicular 

diffusion is ^ lo\ hence we can ignore perpendicular diffusion. For 

q 

low temperatures, ^ 10 K, the mean free path is so short that both 
parallel and perpendicular diffusion are negligible. 

The coefficient of viscosity, , is given by Spitzer (I 962 ) 

T|| = 10 ^ T "^ 2 gm cm sec ^ . ( 3 >l 6 ) 


For our values of density, temperature, and size scale, the Reynolds 
number for the plasma is. 


R 

u 




^ 10 3 


(3-17) 


This large value for R^ means that the viscosity results in a negligible 
dissipation of momentum and energy, i.e. the plasma is approximately 
mviscid In the numerical calculations we introduce an artificial 
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viscosity m order to take care of the possible formation of shocks, 
however we ignore the ordinary plasma viscosity. 

3*5 Magnetic Field Structure 

—4 

In our model it is assumed that all quantities vary along B only, 

hence, one-dimensional equations can be usqd. Justification for this 

assumption is provided by the following arguments. 

The magnetic field in a post-flare region is observed to have 

sufficient strength to dominate the plasma. This may not be true 

immediately after a flare, since, if the plasma thermal energy is 

derived from magnetic energy, then probably they are initially of 

v 

comparable magnitude. However, as the plasma cools from initial 

7 7 

temperatures of ~ 3 X 10 K to < 10 K the magnetic field begins to 
dominate. For our temperatures and densities the field strengths 
reguired are, B ^ 10 G, a value well within the magnitudes observed 
(Section 2.2.1). 

Since the plasma is fully ionized, the conductivity is high 
enough to insure that the "frozen in" condition is valid, i.e. the 
magnetic flux through any area moving with the plasma remains constant. 
In our case the plasma pressure is not sufficient to transport the 
field; hence, the plasma is constrained to move only along field lines 
like beads on a rigid wire. It was shown in the previous section 
that the only significant components of the stress tensor 11 and the 
heat flux q are also parallel to B. Therefore, the equations of 
motion ( 3 * 9 ) reduce to a single equation for mass flux along B, and 
the heat equation (3*10) involves heat flux only along B as well. 
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In order to obtain a one-dimensional model we must assume that 

the relevant variables T, v and n do not vary over some area. A, 

1 ^ 

perpendicular to B. The observed morphology of loop prominences 
(Bruzek, 1964a) suggests that they form along field lines, the loop 
consisting of a flux tube with cross-sectional radius r^ small 
compared to the loop length. For r^ sufficiently large then the 
arguments above imply that the tube is effectively insulated from the 
rest of the corona, since during the lifetime of a loop very little 
heat or mass is exchanged between the loop and its surroundings. This 
will be true when the cooling time t of a loop is much less than the 
time scale for significant transverse conduction, T . The determining 
factor is t/t , where 


T 


1 



K T 



B 2 T 1/2 


-1 

n 



(3.18) 


x 

o 7 

and t s=s 10 sec from observations. Hence for B rj30 G, T pa 10 K, 

and n « 10^ cm ~ = 1 for = 10^ cm. If the tube radius is less 

than 10 u cm. then conduction across the tube is dominant; and, the 

temperature will be approximately constant over the cross-sectional 

area of the tube. However, this means that a significant amount of 

heat may be exchanged between the tube and the surrounding corona. 

Conversely, if the radius r^ is greater than 10^ cm. , then the tube is 

effectively insulated from the outside but significant temperature 

variations can occur across the tube. 

g 

Loop prominences have radii, r^ ~ 10 cm., and so we expect the 
latter case to hold. Hence, there could be small scale irregularities 
across loops, but structures of this size scale, ^ 10^ cm., would be 
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very difficult to observe. In our model we assume that there are no 
small scale irregularities, therefore, the loop is both insulated and 
constant on a cross-sectional area. This assumption allows us to use 
one-dimensional equations. Of course, there is a boundary between the 
loop and the surroundings over which large temperature differences 
occur, however, the discussion above indicates that the width of this 
boundary is small compared to the width of loop prominences Since 
we do not specify the magnitude of the area of the tube m our equations, 
the calculations also apply to the case where the tube area is small 
enough to insure that no small scale irregularities are present, but 
the tube is also insulated from the surrounding plasma. This situation 
might be possible if the field is non-potential. 

Figure 3*1 illustrates the geometry of our model loop. We assume 
a dipole source imbedded a distance D below the top of the chromosphere, 
and consider a flux tube of height above the source R that lies in a 
plane perpendicular to the solar surface. We could easily consider 
the tube to be inclined at some angle from the vertical. This would 
only reduce the effective gravitational acceleration g by some constant 
factor, cos <p, and would not appreciably affect our results. In general, 
loop prominences are not observed to occur at large angles from the 
vertica 1. 

Letting s be distance along the loop measured from the top, we have 
that all variables, including the cross-sectional area A, are functions 
of s only. The effect of the magnetic field on the plasma is accounted 
for by including A and by permitting only heat and mass transfer along 
s. Each flux tube is uniquely characterized by two parameters the 
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size R and the compression factor T which is defined as the ratio of 
field strength at the chromosphere (s = ^ s^) to that at the top (s = 0) 

T = B(s b )/B(s = 0) = A(s = 0)/A(s = s b ) (3*19) 

3*6 One-Dimensional Equations 

Modifying the MUD equations (3 8) - (3*10) to fit our one-dimen- 
sional model and neglecting viscosity, we obtain the following set of 
equations 

f + 5Ti7 h W ‘ 0 ■ (3 - 20) 


n 



V + ^ 
% 



(s) = 0 


(3.21) 


and 


+ v 


3s 


+ E !_(*.) l|_(«|I) + n 2 A(T) . 0 

A 9s A os 3s v ' 


( 3 . 22 ) 


The derivation of equations (3.20) - (3-23), including the viscosity 

term is given in Appendix A. In these equations v is the velocity 

along the field, and g is the component of the gravitational field 
— * 

along B. The variable n is the number density of ions (mostly protons 
in our case), hence the mass density, 


P = (m e + m ) n = n^n , 

the pressure, 


(3.23) 


p = (n + n ) kT = 2nkT , 
r x e p 7 

and the internal energy per ion. 


(3.24) 
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e 


(3.25) 


" I A ’ 3kT 


The area of the tube is given, by A(s), however-, from t-he equations 
above we note that the magnitude of A is unimportant, only relative 
values enter. Therefore, without loss of generality we set the area 
at the top to unity. 


A(s - 0) = 1 (3 26) 

Equations (3-20) - (3*25) are sufficient to determine our four 

variables v, t, n and p as fundtions of t and s. A and g are known 

II 

functions of s, they are completely determined by our choice of 
magnetic field. In Appendix A the dependence of A and g^ on s is 
derived. Defining, 


V = sin 9 , 


where 0 is the angle from the vertical (Figure-3.1 ) . We obtain 


s = | (WSTTT +y== j£n (JT V* + ^3V 2 + 1)} , 


A(V) 


= (1 - v 2 ) 3 /v / 3v 2 + 1 , 


(3.27) 


(3.28) 

(3-29) 


and 


8 1| - 8 s 3v / 


1 - V 
2 

3v + 1 


(3-30) 


Although g^ and A are not expressed explicitly in terms of s, their 
numerical value for a particular s can easily be obtained using 
equations (3-28) - (3»30). Since the heights of the loops we consider 
are small compared to a solar radius, equation (3 -30 ) does not include 
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the effect of a changing gravitational potential with height. 

The only other variables in our equations, the conductivity 
K and the radiative loss rate A, are functions of T only. K is 
given by equation (3*7) and A has been tabulated by Cox and 
Tucker ( 1969 ), Figure 1.1. 
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4 . numerical methods 


4.1 Introduction 

We wish to solve equations (3.20) - (3.22) along with suitable 
initial and boundary conditions. In order to facilitate the computa- 
tion, we first convert to Lagrangean coordinates. There are two reasons 

for doing this- first, we eliminate all "convective" terms (of the 

Sf 

form v second, we are able to treat parts of the plasma of varying 

density with the same accuracy. If thermal instability occurs, then 
xte expect to obtain variations in density of several orders of magnitude 
over the length of a loop. This variation is impossible to follow with 
a grid of fixed intervals m s unless the grid is very fine. With a 
reasonable grid (~ 64 points over the length of the loop), a condensation 
may occupy less than one grid interval. One way to circumvent this 
difficulty is to choose a grid so that equal intervals contain equal 
amounts of p„lasma; the spatial positions of these grid points are no 
longer fixed and will have to be calculated along with the other variables. 

Hence, we change independent variables from s, the distance along 
the loop measured from the top to x(s), the amount of plasma m the 
loop from the top to position s at t = 0, i.e. 
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volume of the loop between the top and the point s. 

A new variable, S(x,t), must now be introduced. It denotes the 

Eulenan position of that slab of plasma labeled "x". Converting 

equations (3.20) - (3-22) to Lagrangean formulation and rearranging, 

/ 

we obtain 


3_ v(x,t) 
dt 


1 <L_ P(x,t) 

A(S) 3x 



(4.3) 


If s(Xlt) - »<*.*) . 


(4.4) 


<L T (x , t ) _ _ p(x,t) d_ (a(S) v(x,t)) _ n (x, t ) A(T) 
dt 3k 3x ' 3k 


2 . 10~ 6 3 _ 

21 • k 3x 


(* 2 - h' 17/2 ) • 


(4.5) 


with 


P = 2knT , (4.6) 

and 

n(x,t) = J^A(S) ||j . (4.7) 

Equations (4.3) * (4.7) are m the form that we use for numerical 
calculations. We assume a grid of equal intervals m x, which is 
equivalent to dividing up the loop into slabs of equal mass. The 
behavior of the loop plasma is then simulated by following the evolution 
of these slabs as determined by equations (4.3) - (4-5). 
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poo a 

4.2 Boundary and Initial Conditions 

In order to determine the correct boundary conditions for the 
system (4.3) - (4.7), we first note that there are actually only two 
distinct dependent variables the position of the plasma S and the 
temperature or internal energy T. Equations (4.3), (4 6) and (4.7) 
can be used to eliminate v, p and n, resulting in one partial differential 

I 

equation of second order in time for S and one of first order in time for 

T. The spatial terms, however, are highly nonlinear. The highest x 

derivative for both S and T are of second order. We essentially have 

a hyperbolic equation for S and a parabolic equation for T m the two 

independent variables x and t. Suitable boundary conditions for such 

a system are Cauchy conditions for S and Dirichlet for T on an open 

region m the x,t plane. In particular, as initial conditions we 
dS 

specify T, S and ^ on t = 0, and as spatial boundary conditions we 
specify T and S at both ends of the loop for all t. 

From the geometry of the loop. Figure 3*1 and from observations, 
it appears reasonable to assume that the loop is symmetric about the 
vertical axis. This implies that T must be a symmetric function of 
x and S antisymmetric; hence we need solve the equations only for 

^ x ^ 0. At x = 0, i.e. the top of the loop, the boundary conditions 
must be ^ = o and S(0,t) = 0. 

The conditions at the base of the loop are determined from physical 
consideration. We are interested m the cooling phase of flares, 
specifically when radiative losses dominate the cooling. During this 
phase we do not expect to have significant heating or evaporation of 
chromospheric material; therefore, we take the temperature at the base 


reproducibility * 
original pass is 
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k 

of our model to be that in the upper chromosphere, T (S^) « 3 X 10 K, 
and to remain constant during our calculation. We do include the 
conduction of heat from the loop to the base; however, we assume that 
this energy can be radiated away by the base without appreciably 
affecting its temperature. 

This situation may appear inconsistent since on the one hand, we 
have that the cooling rate of the base is high enough to radiate all 
of the energy conducted down from the loop yet, on the other hand, the 
rate is so slow that the base temperature does not change during the 
cooling time of the loop. The resolution to this inconsistency is 
m the behavior of the radiative losses at low temperatures. As 
discussed m Section 3*3 , the loss rate decreases very rapidly for 
low temperatures. Assuming the Cox and Tucker loss rate, we obtain 
for n ~ 10^ cm ^ and T ~ 3 X 10 4 K 

k _o _ i 

n e nA(T) ~ 10 ergs cm J s 

This radiative rate is high enough to dissipate an incoming energy 
9 -2 -1 

flux of 10 ergs cm s m a thickness of only 1 km. But at 

T ~ 10 4 K, n g nA(T) ^ 10 ^ ergs cm ^ s \ which implies a cooling time, 
3nkT 5 

T c ~ n~nA ^ 10^ s * Therefore, the base plasma is "trapped" between 

4 n&a 4 

10 - 3 X 10 K; the cooling rate at the lower temperature is too low 

to allow a significant decrease, while the losses at the upper value 

are too high to permit an increase m temperature. The exact profile 

of the base temperature is unimportant to our model, hence, for 

4 

convenience we take it to be constant at 3 X 10 K. 
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As a boundary condition on S, the position of the base is assumed 

to be fixed during the cooling process. This is a good approximation 

since the gravitational scale height of the chromosphere is very small 

4 kT i 

compared to the size of typical loops. At T ~ 10 K, H = ~ 400 km; 

therefore, the chromospheric material will move very little m response 

to pressure changes in the loop plasma. 

The initial conditions for our model loop are somewhat arbitrary. 

Order of magnitude values of the temperature and density are set by 

7 11 -3 

observations, T ~ 10' K and n ~ 10 cm , however, the detailed profiles 
are not known. As conditions on the velocity and pressure, we assume 
that the plasma is initially at rest, 

v(s, 0) = 0 , (4.8) 

and in hydrostatic equilibrium, 

= m^. n g . (4.9) 

Neither of these assumptions is probably valid immediately after the 
loop fills with evaporated chromospheric material, however, we expect 

these conditions to be achieved m a time short compared to the cooling 

3 9 

time, ^ 10 secs. For loop heights of ~ 3 X 10^ cm. the travel time 

2 

across the loop at the sound speed is only ~ 10 s; hence, plasma motion 
will cancel unbalanced pressure gradients before appreciable cooling 
can take place. This motion will m turn be damped by thermal conduction; 

Q 

the dissipation time for wavelengths S 3 X 10° cm. is ;£ 30 seconds. 

g 

Oscillations of long wavelengths, ^ 10 cm., corresponding to a funda- 
mental mode of the tube could persist; however, our calculations 
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(Section 5.4.3) indicate that they do not significantly affect the 
formation of condensations. 

The important parameters of the initial temperature profile are 
the average temperature, which is determined from soft x-ray observations, 
and the heat flux to the chromosphere. Only an upper bound can be 
placed on the latter, the downward heat flux should not be more than the 
total radiative rate during maximum flare brightness. (This argument, 
of course, assumes that a significant portion of the flux is not dissi- 
pated by evaporation.) Since flares are not usually seen m white light, 
we can place an upper bound of ~ 10^ ergs cm ^ s ^ on the heat flux 
into the chromosphere, however, our calculations (Section 5*2), indicate 
that the flux is usually much less than this value. 

Assuming that the heat flux is of the same order as the chromospheric 
7 -2 -1 

HO? flux, F ~ 10 ergs cm s , we obtain for the temperature scale height; 


K 


T 


/l dT 

It ds 


-1 



7/2 


(4.io) 


11 7 

hence = 3 y 10 cm. at T = 10 K. Since this scale height is much 
larger than the size of a loop, only a very slight temperature difference 
over the loop length is needed to carry the flux. Also, the scale height 
has a strong dependence on T. For low temperatures, T ^ 10^ K, equation 

g 

(4.10) yields ^10 cm , which is much less than typical loop lengths. 

From the arguments above, we expect that most of the loop volume will 
7 

be at ~ 10 K, with only a small fraction of the plasma at lower temper- 
ature. Therefore, we approximate the profile by assuming that initially 
the temperature is uniform at ~ 10' K, i.e. 


T (s, 0) = const. 


(4.11) 
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Equation (4.9) can now be used to determine the pressure profile 


p(s, 0) = p Q exp 


f s f 
2kT J g u 1 

0 


s) ds 


(4.12) 


However, for T ~ 10 ( K, the gravitational scale height 


H = > 10 11 cm. 

8 "H 8 !! 


(4.13) 


which is much larger than the loop length. Hence, equation (4.12) can 
be well approximated by 


p(s, 0) = p Q = const. 


(4.14) 


Equation (4.11) and (4.14) imply that 


n(s, 0) = const. 


(4.15) 


Using the initial density profile (4.15) and equation (4.1), we 
can express the initial conditions (4.8), (4.11) and (4.15) in terms of 
x. In this case, the procedure is trivial; all three initial profiles 
are uniform in x. They have the obvious advantage of being very simple. 
Perturbation can easily be introduced into these profiles. There is 

the difficulty, however, of matching them with the boundary conditions 

4 7 

at the base. Since the base is at 3 X 10 K while the loop is at 10 , 

we, essentially, have replaced the steep gradient m T by a discontinuity. 

Therefore, we must modify the boundary conditions at the base. 

The important quantity is the heat flux out the loop. We use the 

following method to determine it. Assume that the loop has been divided 

into a number of equal mass intervals. The position and the temperature 

T^ of the base are known and assumed to be constant m time; also, the 
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mean position and temperature of the plasma in the lowest slab, 

$£ and Tf, are known at some time t. The outward heat flux is then 
calculated by assuming that the flux between and is independent of 
x, i.e. we assume that for fixed t the heat flux across the final slab 
is a constant heat flux determined by its temperature and width. Thus, 


F 


b 




(4.16) 


where we have used the Spitzer conductivity. If the radiative losses 
of the plasma between and are small, then the expression (4.16) 

is valid. If the losses are large, then we are overestimating the 
outward heat flux since a significant portion of it will be radiated 
away before reaching the base. In our calculations, - S^, is usually 
small enough so that the first situation is realized. In any case, the 
maximum error is only of the order of a factor of 2 m the total energy 
losses of the last slab. Also, our loss rate function A(T) is probably 
only of that order of accuracy. Since we are primarily interested m 
the behavior of the plasma at the top of a loop, -where condensations 
form, we do not expect that the use of equation (4.l6) introduces a 
significant error. 

In summary, we have that equations (4.3) “ (4-7) are to be solved 
on the open region, , 


0 £ x £ x , and t ^ 0 , 

m 

with given initial conditions 

T(x, 0) , S(x, 0) and S(x ’ 0) (= v(x, 0)) , (4.18) 
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and spatial boundary conditions 

3_T(0, t) = 0 , T(x , t) a 3 X 10 ^ K , 

3x m 

S(0, t) = 0 and S(x , t) = = const. 

Note that equation (4.4) and the conditions above imply that 


(4.19) 


v(0, t) = v(x m , t) = 0 . 


(4.20) 


4.3 Shocks 

A sticky problem m the numerical solution of hydrodynamic equations 
such as (4.3) ~ (4.5) is the handling of shocks; no fast method for 
treating them exactly exists. We use the standard procedure (e. g. 
Goldsmith, 1970) of introducing the van Neumann-Richtmeyer artificial 
viscosity term m the equations. This term is a pseudo-viscous pressure. 
In equations (4.3) and (4.5) the pressure p is replaced by (p + q) where 
q is defined as 



q 







(4.21) 


where Ax is the grid spacing and a is a dimensionless constant, a ~ 1.5 ' 

2, which determines the thickness of the shock. 

The effect of introducing this term m the fluid equations has been 

investigated m detail (Richtmeyer and Morton, 1967)- A true shock 

2 

shows up as a jump in the values of n, T and v with a width of ~ a 

intervals and with approximately correct values for the size of the 

/ 3v 

jump and the speed of the shock through the fluid. Because —I is 
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small outside of the shock, q does not affect the behavior of the plasma 
away from the shock. In addition, since q is defined to vanish for 
positive velocity gradients, it has no effect on free expansion. Tests 
on the artificial viscosity by various authors (Richtmeyer and Morton, 
I967) indicate that, m general, it seems a fairly accurate (and fast ) 
way of accounting for shocks. 

There are two places m our model where we expect that shocks might 
form, and hence, q becomes important. One is at the boundary of a 
condensation where large pressure gradients form. The size scale h 
of a knot is determined by a balance between radiative cooling and 
conductive heating, 

h w 10~ 3 T 7/ ^ . (4.22) 

nsA(T) 

The time scale for growth of a knot is the cooling time. 



(4.23) 


In order to maintain pressure equilibrium, plasma must move fast enough 
to keep up with pressure changes, 


v « h/r « 4 X 10 12 v/a""t^ . 

C 


(4.24) 


As long as these velocities are damped quickly (i.e. in a time short 
compared to the cooling time) then they remain subsonic. 


M = ~ = 4 X 10 8 ,/A t 1 ^ , 


(4.25) 


4 6 

and M < 1 for 3 X 10 £ T s: 5 X 10 K. However, since conduction is the only 
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dissipative process present, then these velocities will not be damped because 
whenever condensations occur, the radiative cooling time must be short 
compared to the conductive cooling time. Hence, velocities generated 
while the plasma is hot will persist and result m shocks. At T ~ 10^ K, 
equation (4.24) gives v = 30 km s \ which is supersonic at T = 3 X 10^ K. 

In Section 5*4.2, we examine a simple analytic model for the growth of 

condensations and find that shocks may indeed form at the boundary 

« 

between the cold and hot regions. Our numerical calculations also 
support this result. 

Shocks also form when cool, freely-falling material impacts onto 

7 -1 

the chromosphere. Downward velocities on the order of ~ 10 cm s 

/ 

are observed; these are well above the sound speed at chromospheric 
temperatures. In our model, the artificial viscosity q dissipates the 
streaming energy of the falling material. The usual dissipation mech- 
anisms, conduction and viscosity vary as high powers of T and are 
negligible at low temperatures. However, the velocities of, the cool 
plasma are obviously damped by some means since this plasma is not 
observed to bounce back into the corona after impacting the chromosphere. 

We believe that the most likely source of this dissipation is shocks. 

4.4 Computer Code 

In order to obtain numerical solutions for the system (4.3) - (4.7), 
including the artificial viscosity, we use a two-step, second order 
difference scheme similar to the ones described by Richtmeyer and 
Morton, 1967* We have three partial differential equations to solve 
for the variables v, S and T. They can be written as: 
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■where V is a three component vector (v, s, T) and F is a nonlinear 
function of V and its spatial derivatives. The pressure and number 
density are given as explicit functions of V by equations (4.6) and 
(4,7); hence, they can be obtained directly. Assuming that V(x, t) 
is known, V(x, t + At) can be obtained by a Taylor series expansion. 


V 


= V 

t +.At t 


= Vt + At IF i + 


t 


At 

2 


cJF 

at 


+ 0 (A ) 


(4.27) 


But, 

F |t + At= F |t + 4 |lf|t + 0 ( it2 ) 

2 

hence 


(4.28) 


V | t+ At- V |t +4tF |t + 4 JL + °( 4t3 > ' 

2 

From equation (4.28) we note that ^ | t + LS needed only to first 

2 

order accuracy in At; therefore we first generate provisional values 

^ * 

of V i . defined by: 

|t + At 7 

2 

A £ O 

V |t+At SV |t + - F |t-' , |c + At + °( 4t > • fr-30) 

2 2 

and define 

F |t + A| SF ( 7 |t + A|j= F ( V | t + A|j + °( 4t2 ) ' 


Since F is a smooth, well-behaved function of V, we can expand the 
right-hand side of equation (4.31) to obtain 



(4.32) 





4~ A t 
2 


t + 


At 

2 


+ 0(At ) 


The scheme consists of first generating V| fc + from (4.30) and then 

~ 2 

calculating F j t + then V j t + is obtained from 

~ 2 

v |t + 4t = v |t + “*!<: + At + 0 < At3 > • ( 4 - 33 > 

2 


Equation (4.33) follows immediately from equations (4.29) and (4.32). 

The boundary conditions S(x, 0), v(x, 0) and T(x, 0) are used to 
initiate the calculation and the spatial boundary conditions are used 
to compute the velocities and heat flux at the ends of the space grid. 
To facilitate the computation, we rewrite equation (4.7) as 


n(x, t) 


3_AI (s ) 
S x 


-1 


(4.34) 


where 

Al(s ) = f A(s ) ds 

4 ) 

Three known functions of S, g^(s), A(s) and Al(S), appear in the 

1/2 

equations and two of T, A(T) and T . Rather than evaluate these 
five functions directly (this requires too much computer time), we 
first tabulate them for ~ 800 semi- logarithmically spaced values of 
S and T, and use linear interpolation to estimate their value for a 
particular S or T. The estimated values are accurate to better than 


(4.35) 


The fineness of the x grid. Ax, and the time increments At are 
determined by the requirements that they be smaller than length scales 
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and time scales of interest and that the difference scheme remain stable. 

The equations that we wish to solve are highly nonlinear. Unfortunately, 
very little is known about the stability or convergence of finite difference 
schemes for such systems. The only accurate method known to determine 
stability is, of course, trial and error. However, we calculate stability 
criteria for several linear approximations to our system and hope that 
these restrictions are sufficient to assure convergence. 

If we neglect the heat loss terms, the area factor and gravity, 
our equations reduce to those of a perfect fluid. Our scheme then 
becomes equivalent to the Lax-Wendroff scheme (Richtmeyer and Morton, 

1967)j which has as its stability criterion 


AS 


At, < Mm 


( jv | + C ) » 

J J 


( 4 . 36 ) 


where AS^ is the Eulenan width of the grid spacing Ax at the jth grid 
point, and is the speed of sound at this point. The index j runs 
over all the intervals m the x grid. If we approximate the heat 
equation (4.5) by a simple diffusion equation 


dT iA 3 2 T 


( 4 . 37 ) 


then the stability criterion is 


(A*r 

Atg < Mm — £- L 
J J 


( 4 . 38 ) 


Since the radiative term is also important, we introduce another restriction 


At < Min p h , 
J J J 


( 4 . 39 ) 


where p^ is the pressure at the jth grid point and £ is the radiative 
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losses there. Combining these three criteria, we adjust the time 
increment at each time t so that 


At < Min (ht v At 2 , At 3 ) (4.40) 

The spatial grid spacing. Ax, is kept fixed throughout a run in 
order to assure proper centering of the difference equations. From 
symmetry arguments. Section 4.2, only half the loop need be included. 

Q 

We use 33 points for the x-grid, hence, for a loop of height 2 X 10 cm,, 
the average width of each grid interval is 700 'km. The width of a 
particular interval, i.e. of a particular slab of plasma, depends on 
the density and area at that point; from equation (4-7) we have 


As 


Ax 

nA(s ) 


(4.4l) 


Equation (4.4l) implies that m regions where the density is high or 
the area large, such as in a condensation at the top of the loop, the 
Eulerian spacing of the grid is finer. As mentioned previously, this is 
the mam advantage to using a Lagrangean coordinate system. 

For the temperatures, densities and sizes of typical loops, and 
with 32 slabs for the x grid, we find that the time increment At, 
as determined by equation (4.41), ranges from At .05 - *25 sec. To 

7 

follow a typical loop from the 10 ,K temperature level until the time 
when most of the plasma has cooled and fallen onto the chromosphere 
requires 1 - 3 minutes of computer time on Stanford's IBM 370/168. 

A large fraction of this time is taken up during the initial conductive 
cooling stage of the loop when criterion (4.38) dominates. Unfortunately, 
the time for a run, t r , varies as a high power of the grid spacing. If 
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equation ( 4 . 38 ) is the important criterion, then t r °= (Ax)~^, thus, we 
are severely limited in the possible fineness of the x grid. Details 
of our difference scheme along with a copy of the code are given m 
Appendix B. 
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5. COOLING OF FLARE LOOPS 

5.1 General Evolution of Loops 

From our computer calculations the following picture emerges of 
the evolution of a post-flare loop. The evolution can be divided into 
three distinct stages characterized by the behavior of the temperature 
profile. In the first stage — immediately after flare heating ends — 
the temperature is high enough so that conduction to the chromosphere 
dominates the energy losses. We expect that evaporation may also be 
an important cooling mechanism during this phase. During the first 
stage, the plasma is stable against the growth of thermal perturbations; 
hence, a smooth temperature profile develops. The scale height is so 
large that the loop is approximately isothermal, and the cooling rate 
varies as a high power of the temperature. These results also hold 
if evaporation is included (Antiochos and Sturrock, 1976^)* 

As the temperature decreases the effectiveness of conduction drops 
rapidly due to the strong temperature dependence of the thermal conduc- 
tivity, equation (3.7) > and eventually radiation begins to dominate the 
energy losses. The loop then enters the second stage of its development, 
during which radiation is the primary cooling mechanism. It is during 
this phase that thermal instabilities can occur; hence, the evolution 
of the loop during the second stage is highly complex involving large 
temperature gradients and plasma velocities. We find that the details 
of the evolution, in particular whether condensations form, are highly 
sensitive to the initial configuration of the loop plasma when it first 
enters the radiative cooling stage. 
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Since the density m post-flare regions is high, radiative losses 
are large and the temperature throughout a loop quickly drops to chromo- 
spheric values , T ~ 10^ K. Once it reaches this value the temperature 
changes very little during the subsequent development of the plasma 
because the heat losses at low temperatures are insignificant. We 

define this phase of the evolution, during which the plasma remains 

4 

isothermal at ~ 10 K, as the third stage. In this phase large density 
gradients and velocities are generated due to pressure forces and the 
effect of gravity. 

In our model, the end point of the evolution is the establishment 
of a gravitational scale height atmosphere at 3 X 10^ K. This does not 
agree with the true solar atmosphere since we have not included coronal 
heating m the calculations. For our large initial densities, the 
energy loss rate of the loop plasma is much larger than the estimated 
value of the coronal heating rate, ^ 10° ergs cm s , (Moore, 1972), 
However, once the loop plasma cools to chromospheric values and begins 
to fall onto the chromosphere, then the density in the upper part of 
the loop decreases until coronal heating once again becomes important. 
Therefore, we terminate our calculations when the density at the top 
of the loop drops by more than an order of magnitude from its initial 
value . 

The time scales for the 'three stages of evolution are given by: 

(i) the conductive cooling time at high temperatures, T ~ 10^ K; (n) 

g 

the radiative cooling time at T ~ 10 K; and (m) the free-fall time 
along the loop. For typical loops, n = 5 X 10 10 cm"^ and H= 3 X 10^ cm, 
these times are 
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* 


(5.1) 






T = 2/^-2- « 10 2 s 
T ii 2. ~ iU S » 

n A 


(5.2) 


an) t 1±1 = A /-~-«io 3 s . 


( 5 . 3 ) 


We note that the time scale for stage (li) is shorter than the lifetimes 

of condensations in loop prominences ~ 10 minutes; hence, the bright HQ' 

knots cannot be due simply to a temperature difference along the loop, 

instead they must represent a density difference. Our calculations verify 

this hypothesis. We find that due to the thermal instability of the 

plasma, large temperature differences form; but, these die out quickly 
2 

^10 s. However, the temperature gradients give rise to density 

q 

enhancements which persist for long times ~ 10 s, because the speed of 

sound is small at low temperatures. For example, m order for a 

4 g 

condensation at 10 K to expand to a size of 10 ^ cm at the sound speed, 
the time required is 10 min., which is of the same order as the life- 
time of Hq knots. 

In summary, we find that the evolution of post-flare loops is in 
three stages* 

(i) conductive cooling (with or without evaporation) at high temperatures, 
~ 10 7 K, 

(ii) radiative cooling at lower temperatures, ~ 10^ K, 

(m) the dynamics of isothermal plasma at chromospheric temperatures, 
li 

~ 3 x 10 K. 


The stability of a particular loop is, of course, determined by the 
* 

second stage. In the rest of this chapter we describe each of these 
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stages in detail, and determine under what conditions will thermal 
instability significantly affect the evolution of loop plasma. 


5.2 Conductive Cooling Stage 
5.2.1 Conductive Damping Rate 

During the initial phase of flare cooling the temperature is high 
enough so that the conduction term dominates in the dispersion relation, 
equation (l.l4). Therefore, we expect perturbations to be damped on a 
time scale. 


T 


D 


*T(f) 


(5.4) 


since t-q is the time scale for temperature differences over lengths 

X/2 to be significantly reduced by conduction. We tested this result by 

examining the behavior of a perturbation of small wavelength but with 

a finite amplitude. The perturbation was chosen to be initially 

isobaric with an amplitude of 9*75$ and a wavelength \ = 2 X 10^ cm. 

The loop was assumed to be due to a point dipole source, with a height 

above the chromosphere of 3 X 10^ cm and a compression factor T = 10. 

For this choice of height and T, the loop length turns out to equal 

4 wavelengths of the perturbation. The initial average temperature 

and density assumed were T = 10^ K. and n = 5 X 10^ cm ^ . 

In Figure 5*1> the evolutions of the amplitudes of the temperature 

and density variations are shown. We note that the damping time for 

the temperature perturbation is ~ 5 sec . , which is in good agreement 

with eq. (5.4); however, the evolution of the density is very different 

from that of the temperature. The initial response of the density 

* 

profile lags behind that of the temperature due to the inertia of the 
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Figure 5.1 EVOLUTION OF THE AMPLITUDE OF A 
TEMPERATURE AND DENSITY PERTURBATION THAT IS 
DAMPED BY CONDUCTION. The solid line refers 
to the density perturbation and the broken 
line refers to the temperature perturbation. 
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plasma. The maximum velocity generated by the resulting pressure 

' 6 -1 

gradients is ~ 3 X 10 cm s , which is an order of magnitude less 
than the sound speed. 

By t = lh sec. we note from Figure 5*1 that the amplitude of 
the perturbation has dropped by over 60 %, whereas, the thermal energy 
of the loop has decreased due to conduction out the base and radiation 
by less than V } 0 . The short damping time implies that we are justified 
in assuming an isothermal and static initial state for the loop plasma. 
All small scale gradients due to the initial evaporation of chromospheric 
material by flare heating decay quickly compared to the cooling time. 

It is possible that large scale motions, on the order of the length of 
the loop, and of small amplitude persist for time scales on the order 
of the cooling time. However, our computer calculations indicate 
that velocities which are small compared to the sound speed have 
negligible effect on the evolution and stability of flare loops. 

If evaporation cooling is important, then large plasma velocities 
may be driven by this process. We can estimate the magnitude of these 
velocities — the downward heat flux must be approximately equal to the 
mechanical flux upwards, 


pv « F . ( 5 . 5 ) 

8 Q -2-1 

Taking the heat flux to be of order 10 - Kr ergs cm s for large 

12 -3 

and small loops, respectively, and the pressure ~ 10 - 10 ergs cm , 

we obtain, 

7 -1 

v w 10 cm s . ( 5 . 6 ) 
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Since the sound speed at 10 K is 50 km s , the Mach number is 
M = v/c = .2, and the ratio of kinetic to thermal energy is only 
.04. In addition, we expect that since the heat flux F « the 

Mach number M “ T^, if the pressure is constant, or M “ T^, if density 
is constant. In either case we note that the Mach number varies as a 
high power of T; hence, even if the evaporation velocities are 
comparable to the sound speed when the loop temperatures are high, 
these velocities will be negligible when the temperature drops low 
enough so that radiative cooling dominates. 

5.2.2 Static and Evaporative Cooling Rates 

The first stage cooling of flares has been investigated by several 
authors Culhane et al. (1970), Strauss and Papagiannis (1971), and 
Zaumen and Acton (197^- )• la their models they assumed that a loop 
has uniform cross-section; however, Antiochos and Sturrock (1976a) 
point out that for magnetic fields typical of flare loops (Rust and 
Bar, 1973), the variation of area along the loop has a significant 
effect on the cooling. Antiochos and Sturrock have calculated the 
evolution of the temperature and density profiles both for static and 
evaporative ( 1976 b) models; they assume that gravitational and radiative 
effects are negligible and that all velocities are small so that 
pressure gradients are unimportant. With these assumptions the MHD 
equations simplify sufficiently so that analytic solutions may be 
obtained. 

For completeness we state the important results of the calculations 
by Antiochos and Sturrock below. In the case of static cooling of a 
dipole loop, the temperature profile is given by 
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-ojr r / ^ci - <)\ ~i 

T(e,t) = T(o,o)(l + t/r) 2/5 1 -( ~) , ( 5 . 7 ) 

L V\d - v 2 )/ J 

where V has been defined, in equation ( 3 . 27 ) and the subscript "b" refers 
to values at the base of the loop. The cooling time t is given by 

V 2 

T = 10 6 p(o) (t(0,0))~ 7/2 R 2 , (5.8) 

(l - s>lf 


where R is the radius of the loop. Figure (3.1). The parameter v b is 
directly related to the compression factor T — the dependence of t on 
T is given in Figure 2 of Antiochos and Sturrock ( 1976 a). For the 
case of evaporative cooling, the results are similar except that the 
temperature varies as (1 + t/T) 2/ ^ 7 and at low temperatures the 
heat flux varies linearly with T rather than being constant as m the 
static case. 

o 

We ran our code for a dipole loop with T = 10 and R = 2.5 X Kr cm, 

m order to compare our model with the results above. The initial 

7 10 -R 

temperature and density were chosen as 1.05 X 10 K and 2.5 X 10 cm , 

respectively. With these values the cooling time given by equation 

/ w-U"'5/2 


(5-8) is 200 sec. In Figure 5«2 we plot 


\T(0 )} 


- 1 versus t. 


where T is the average temperature of the loop in our computer run. 

We note that the graph is linear in agreement with equation (5»7) and 
that the cooling time is I 70 sec. which also agrees with theory. The 


loop size and plasma density were deliberately chosen to be small so 
that radiative cooling would be negligible. We find that the plasma 
loses 32 $ of its thermal energy m 270 sec., with 88$ of the loss due 
to conduction and the rest via radiation. 
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Figure 5.2 COMPARISON OF ANALYTIC CONDUCTIVE 
COOLING MODEL WITH OUR NUMERICAL RESULTS 
PLOT OF (T/T )-5/2-l VERSUS t. 
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These results indicate that in the limit of small radiative losses 

and plasma velocities our numerical code does reproduce the analytic 

calculations. Of course, this cannot be used as an argument to prefer 

static over evaporative cooling since our model assumes no evaporation. 

In fact, we expect that evaporation is probably the dominant cooling 

process at high temperatures because the heat flux predicted by both 

our code and the analytic model is high. The initial heat flux is 
9 -2 “1 

2 X 10 ergs cm s , which is much too large to be ahsorbed by the 
chromosphere. We believe that the initial cooling is accurately 
described by the model of Antiochos and Sturrock (1976b) and that the 
later stages, when radiation is important, are accurately described by 
our numerical code. 

From the arguments in the previous two sections, we reach the 
following conclusions. The evolution of the loop density and, hence, 
the stability of the plasma, are not uniquely determined by the initial 
temperature and density since evaporation may occur. Even if evaporation 
is present, however, the first stage cooling process will bring the 
loop into an approximately isothermal and static state. Given the 
density at lower temperatures, when radiation losses dominate, 
the stability may be determined. In the next section we discuss 
necessary criteria on the density and temperature profiles for the 
formation of condensations. 

5.3 Instability Criteria 

Using our model of flare loops, we have determined criteria under 
which thermal instabilities in a coronal plasma can grow. The term 
instability is not strictly applicable to our model since the plasma 
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is never in an equilibrium state but is continuously evolving, however, 
we use it to denote situations m which the cooling is highly non- 
uniform. By the term "unstable loop", we mean that a condensation 
observable m Ha and of small size scale forms. Specifically, we 
require that the condensation be at least an order of magnitude brighter 
m Ha than the rest of the loop and have a lifetime of the order of 
minutes. These condensations correspond to the bright knots m Ha 
that are often observed to form at the tops of loop prominences 
(Section 2.2.4) . 

We find that the following conditions are necessary for a loop 
to be unstable (according to the definition above)* 

(1) Since both the growth time of thermal perturbations and the lifetime 
of a loop are of the order of the radiative cooling time, only pertur- 
bations of finite amplitude can lead to observable effects. However, 

we find that the amplitude required is quite small, of order 5 $* 

(2) From linear theory (Section 1.2), we require that the growth rate 
due to radiative cooling dominates the damping rate due to conduction. 
Hence from equations (I.I 9 ) and (1.20) we require, 

nX ^ 4 X lo " 3 T 7 ^ A " 1/2 cm -2 , ( 5 . 9 ) 

where X is the wavelength of the perturbation and A is the radiative 
loss function. 

(3) In order for a density perturbation to grow, the cooling time must 
be long compared to the speed of sound travel time across the pertur- 
bation; hence from equation ( 1 . 23 ), 

nX ^ 4 X lO " 11 T 3/2 A " 1 cm " 2 . (5.10) 
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(4) Additionally, we find that for temperatures and size scales relevant 

to loop prominences, the number density m a loop must be greater than 
10 ' 

approximately 10 cm . For densities less than this value the radiative 
cooling time is so long that gravitational effects and conduction to 
the chromosphere dominate the cooling process, stopping the growth of 
condensations. Also, for low densities coronal heating becomes important 
and, hence, our model is no longer accurate. 

The conditions enumerated above agree well with observations that 
loop prominences occur only after flares when the coronal density is 
enhanced and only on large loops of size scales ^ 10^ cm. In Figure 

r 

(5.3) we graph the relations (2) and (3). Inequality (2) is valid 
m the region to the left of the broken line labeled (n) and inequality 
(3) is valid to the right of the line labeled (m). In the region 
between line (n) and (in) isobanc perturbations will grow. If 
conditions (l) and (4) are satisfied as well, then observable conden- 
sations will result. 

5.4 Radiative Cooling Stage 

5.4.1 Comparison with Linear Theory 

We define the second stage of the cooling process to be that phase 
during which radiation is the dominant energy loss mechanism. At the 
onset of instability, we expect that the linear theory may approximately 
describe the evolution. In order to compare our nonlinear model with 
this theory, we followed with our code the evolution of an initial 
perturbation in the temperature and density profile of a loop. The loop 
was chosen in each case to be that due to a dipole field with a compression 
factor, F = 10. The perturbation was chosen in each case to be initially 
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isobanc, with an amplitude of 10$ and a wavelength equal to the length 
of the loop, i.e. we assumed as initial conditions 

T(x,0) = T(0,0) jl - .1 cos(f*)j (5.11) 

and 


p(x,0) = p(0,0) 


( 5 . 12 ) 


where ± ^ are the Lagrangean positions of the base paints of our 
model loop. 

From conditions (5.9) a nd (5.10) 5 the important parameters that 
characterize the stability of a loop are the temperature and the total 
number of electrons, nA, where A is the wavelength of the perturbation. 

In each of the runs the values of T(0,0) and nA were chosen so that 
condition (5«9) was definitely not satisfied and the perturbation was 
initially damped by conduction. However, as the plasma cooled and T 
decreased, (5-9) eventually became satisfied and at a critical temperature 
, the amplitude began to increase. Figure 5-3 shows the value of 
as a function of nA. There is some ambiguity in the exact value of 
since in no case did the temperature profile retain the simple form 
given by equation (5.11). Instead, conductive losses to the chromosphere 
and small but non-negligible velocities altered the profile. Therefore, 
we assumed that the temperature and density variations were due to two 
parts , 


AT 

T 


An 

n 



(5.13) 
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where ) and ( — ) indicate variations due to adiabatic motions of the 
\ T/a \ n/a 

plasma and indicate the amplitude of the isobaric pertur- 

bation. Since for adiabatic motions 


AT 

T 


)a'f(^)a • 


(5.14) 


the perturbation amplitude is given by 


/AT\ / An\ 3 AT 2 An 

( tJi ~ " ( nji 5 T " 5 n 


(5.15) 


Hence, we define T c as the average temperature in the loop at which the value 
of , as given by equation (5.15) 5 began to increase. 

In general the results of our computer calculations agree well with 
theory. We find that T c is approximately independent of n or \ as long 
as the product nl is constant. From Figure 5.3 we note that the effect 
of the nonlmearities m our model is to increase the stability of the 
loop especially at high temperatures. This is to be expected since for 
the initial perturbation with the longest possible wavelength, and, 
hence, the fastest growth rate (equation 5*H)> the temperature maximum 
occurs at the base of the loop. Conduction cooling to the chromosphere 
decreases the temperature near the base, thereby decreasing both the 
amplitude and the length scale of the temperature variations. We find that 
loops with a large compression factor T can support a perturbation with a 
relatively larger wavelength because heat losses to the chromosphere are 
reduced; hence, they are more susceptible to instability, i.e. T in- 

V 

creases with T. Once instability sets in, however, the evolution of the 
plasma is insensitive to T. 

In all cases we find that pressure balance is maintained for only 
a short, initial part of the cooling. As can be seen in Figure 5 - 3 > the 
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temperature T^ occurs near the line labelled at which condition (5*10) 
breaks down and the perturbation can no longer remain isobaric. There- 
fore, significant velocities soon develop. In the next section we 
investigate the growth of these velocities using an analytic model of 
a growing condensation. 

5*4.2 Analytic Model of a Condensation 

We have constructed a simple analytic model to follow the growth 
of a condensation into the nonlinear regime. The model is valid during 
the initial growth when pressure equilibrium is maintained, i.e. 

Vp sa 0 . (5.16) 

Equation (5.16) will be valid as long as the speed of sound and the 
gravitational scale height are large compared to relevant velocities 
and size scales. This model can be used to calculate the temperature 
at which velocities resulting from the difference m cooling rates 
between the cold and the hot plasma become appreciable and, hence, 
pressure equilibrium breaks down. In addition, we can obtain estimates 
of the magnitude of this velocity and determine if shocks form. Of 
course, these results can also be obtained using our full computer 
code, however, we find that the analytic model is useful both to 
check the computer calculation and to clarify the important features 
of the nonlinear instability. 

We assume that the temperatures are low enough so that heat 
conduction is negligible compared to radiative losses. Therefore the 
hot and the cold regions. Figure 5*4, may be taken to be separated by 
a boundary of negligible width. Within each region the temperature 
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Figure 5.4 GEOMETRY OF OUR ANALYTIC 
CONDENSATION MODEL 
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and density are assumed to be approximately uniform and given by, 

^T c (t), n c (t)j and ^T^(t), n^(t)j where the subscripts "c" and "h" 
refer to the cold and hot plasma respectively. Since the two regions 
cool radiatively at different rates, the position of the boundary must 
change in order to maintain pressure balance, equation ( 5 . 14 ); however, 
since conduction is neglected no material moves across this boundary. 

From Figure we note that, 

j& c (t) + ^ n (t) = L/2 = const. , ( 5 . 17 ) 

thus the volume of the condensation is given by 

A (t) 

v c (t) = 2 J c A(s) ds , ( 5 . 18 ) 

0 

and of the hot plasma 

v h (t) = v T - v c (t) 7 (5.19) 

where V T , the total volume of the loop, is constant. Differentiating 
equation (5*18) we obtain the rate of change of the volumes: 

fe - - If V 11 ) - 2A ^<> • (5-2o) 

where v is the velocity of the boundary, 

v - ~ i c (t) . (5.21) 

Since the temperature and density are uniform over each region, we can 
integrate the mass equation over the volume of the condensation 
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X d¥ n c (t) d3y + X n c' (t) v • V d 3 V = 0 . (5.22) 

c c 

Using the divergence theorem on the second term of equation (5.22) and 
substituting equation (5.20) yields, 


v c (t)| J n c Ct) + „ c (t)|jV c (t)=.o , 


(5.23) 


therefore, 


n (t) V c (t) = const. 


(5.24) 


Equation (5.24) is equivalent to the assumption that there is no mass 
transfer between the hot and cold regions. This can only be true if 
the cold plasma cools faster than the hot, i.e. the temperature must 
be in the unstable regime of the cooling curve. Once the temperature 

i 

enters the stable regime, the profile must eventually become uniform 
throughout the loop. 

The heat equation can be integrated m a similar manner as the mass 
equation; yielding 


■3- v ^ + = -An 2 T _a V 

2 c dt 2 ^ dt c 0 0 0 c 5 


(5.25) 


where we have assumed a radiative loss function of the form. 


£(n,T) = A o n 2 T~ a . 


(5.26) 


Defining, 


4 (t) = V (t)/V_ 
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and 


V't) = \(tt/v T , 


(5.27b) 


we obtain the following set of equations 


* c (t) + ‘^(s) = 1 

“c,h (t) - N o,h - COM - 

<L P , i d X (t) _ 

J '-J 11 — 


2 c,h dt 


2 P dt 


A.n 2 , T““ A , 
O c,h c,h c,h 


( 5 . 28 ) 


p = 2kn , T , 
r c ,h c ,h 


The seven independent equations in the system (5.28) are enough to 

uniquely determine the variables (n T X ), (n, T, X, ) and p for 

c ^ c ^ c in 3 in j n. 

given initial conditions. 

The quantity that is of interest is 


5(0 s T c (t)/T h (t) , (5.29) 

the ratio of cold to hot temperature. (The amplitude of the perturbation 
is equal to (l - %)/ (l + §).) After some manipulation, the two heat 
equations m ( 5 . 28 ) can be decoupled and integrated to yield. 


*(S) = P(^)t/T 


(5.30) 


where 


8 (ar+l ) 2tf-l 

P(C)=/[l - ^2j'3(«2) (n+ „3 4t 

r -1 5(0/+ 1) 2&- 1 

pcv-fL 1 - cT 3( ® +2) c 2 <*+v 3 


(5.31) 

( 5 . 32 ) 
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T 


» 


(5.33) 


eIo) 

A 0 n c(°) T J°) 

and 

j 

h=n h /n c . (5.3^) 

There are two adjustable parameters in our solution, Tj, the ratio of the 
total number of cold particles to hot and, the initial temperature 

ratio §(0). We are interested m initial perturbations of small 
amplitude and large wavelengths; hence, 

.9 ^ l Q < 1 and Tl« 1 • (5»35) 

For example, an initial perturbation with an amplitude of 5$ and a 
wavelength equal to the length of the loop corresponds to = .9 and 
T| — 1. 

In deriving equations (5.30 ) and (5*31), we assumed that the 
exponent, ff, of T m the radiative loss function (equation ( 5 . 26 )), is 
a constant. From Cox and Tucker's ( 1969 ) results. Figure 1.1, we 
note that Ot is actually a function of T; however, we can approximate 
the cooling function by segments with constant & as m Moore and Fung 
(1972) and Goldsmith (1971). Since the cooling function is probably 
only accurate to at best a factor of two, the approximation of constant 
01 will not affect the general results of our model. 

The time dependence of ? can be determined from equation (5*30) and 
(5.31); however, they break down for Of = .2 as can be seen from the form 
of (5.31). This is to be expected since we assumed that the loop is 
unstable to the growth' of the condensation, but for a cooling function 
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of the form given by equation (5.26), Field's ( 1965 ) criterion, 


d£ n 
3T T 3n 


(5.36) 


is not satisfied for & £ -2. Thus, our solution is only valid for 
& > -2. For certain & equation (5*31) can be integrated directly 



(5.37) 


From the Cox and Tucker curve we find that 1/2 ^ Ol 1 ^ 2 in the unstable 

5 7 

temperature regime, 10 ^ T ^ 10 . 

In Figure (5.5) we plot § as a function of t/T for each of the 
three values of & in equation ( 5 * 37 ) and for T| = 1 , = .$)• We note 

from Figure (5.5) that for ^ ^ 10 \ £ 13 a very rapidly decreasing 
function of t. Since the plasma tries to maintain pressure balance 
between the hot and cold regions, we will find that the sharp decrease 
in 5 implies equally rapid changes m the density and hence, supersonic 
velocities result. The results in Figure (5*5) are not sensitive to 
the value of T| or Defining t^ as the value of t when § = 0, we 

find from equation ( 5 * 30 ), ( 5 . 32 ) and ( 5 * 37 ) that for GK = 1 / 2 , = 1 , 

independent of T) or For & = 2, t^ -* const, as T] -* 00 or 0, or as 

-> 1. For Qt — -1, t^ diverges for T| -* 0 or ^ -» 1, but only logarith- 
mically. (Also, we are mainly interested only m positive values of a.) 

To obtain the time dependence of the other variables, we express 
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Figure 5.5 EVOLUTION OF THE TEMPERATURE RATIO, F FOR 
VARIOUS Ot. The initial amplitude = .9 and the 
ratio of hot to cold plasma, T\ = 1. 




them in terms of §* 


= §/-m + s ) 


(5.38) 


p/Pq = 


Vl + 


ri + 


i) 


5/3 A - f +2 \ 3 ^ +2) 


1 “ l 


9+2 


0 


T/T_ = 


p K 


C C0 *>0 X cO ’ 


n /n~ = ^ 

c cO cO c 


(5.39) 

(5-40) 

(5.4l) 


The equations for and are identical to (5.4o) and (5.4l) except 
that "c” is replaced by ,r h". 

Our model breaks down when the velocities generated by the unequal 
cooling rates become comparable to the speed of sound. We calculate v, 
the velocity of the boundary between the two regions. Using equations 
(5.20), (5.27) and (5.38) yields, 


L <A> W dg 

2A <+> („+ i ) 2 dt 


inhere we have defined, 


2 

<A> = 7 

J_x 




Defining, 


v 


th 


L_ 
2t ’ 


(5.42) 


(5.43) 


(5.44) 


and using equations (5«30) - ( 5 * 32 ) we obtain. 
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8 a+il 


2a- 1 


2 <A> [l - 


th 


5 A(A C ) 


[1 - 


g 0H-2 1 3(Qf+2) ^ +2 Ti(n+ 3 

5fe+l) 2&+5 

^. 2 ] 3 (» 2 ) ( „ + s) 3 


C5- 1 *?) 


If the loop is of constant cross-section then, of course. 


<A> 

A(A e ) 


= 1 . 


In general, we expect this factor to be less than unity since A increases 
as X (and hence g) decreases. For example, a dipole loop with a 
compression factor F = 30 and a volume for the condensed region, X c = 1 / 10 , 


implies y = .46. Therefore, the variable area factor acts to decrease 

the velocity, improving the validity of our model. A lower limit on 
the size of the area ratio is 


<A> Afs = L/2) _ r -l 


A(A C ) A(s = 0) 


however, this ratio is usually much larger than T 

<A> 

discussion we assume that ' y « 1 . 


-1 


(5.46) 

In the subsequent 


We have plotted in Figure 5-6 — — as a function of 5 f° r 13=1 


v 


th 


and 


2^ = .9, and for several values of 01. Assuming typical loop parameters, 


6 lo -3 Q c 

T ~ 5 * 10 K, n ~ 5 X 10 cm , and L ~ iQr cm implies that ^ 10; 

V th 

hence, our model will begin to break down at approximately this value for 

V /-V 

. From Figure 5.6 we note that = 10 corresponds to a value of 

V th v th 

g ranging from .4 - .2 for .5 ^ & £ 2. These values of g correspond to 

a perturbation amplitude of 43% - 67 %. The temperatures, T and T, , can 

c n. 

be calculated for a particular g using equations (5.33) - (5.40). In 

Figure 5 . 7 , T /T and T./T, - are plotted as functions of g for Ti = 1 and 
c cO h hO 

^ = . 9 . We note that for the values of g at which our model begins to 
break down, the condensation temperature has decreased by approximately 
a factor of 6 from its original value while the hot temperature has 
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dropped by a factor of 2. Since the plasma must cool over several decades 
in temperature (from 10 ' K to 10 K), we expect that large velocities 

(and probably shocks) will be present during most of the cooling. 

It is instructive to examine the dependence of these results on 
the initial amplitude L. From equations (5*38) - (5.40) we obtain. 




5 

3(o+2) 


1 - S 


8+2 





(5.47) 


Letting 8=1-^ and e = 1 - §, equations (5.45) and (5*47) imply that 
as 6 -* 1 , 


and 


e _^\-3(»2)/5 6 



2 <A> T[ 

5 A <V (mf 


(&+ 2 ) 


T C V(8o+ll)/5 

T~ / 
c O' 


6 . 


(5.48) 


Therefore, if the initial perturbation amplitude 6 is small enough, signif- 
icant temperature differences or velocities do not form even if the plasma 
cools considerably, i.e. T c /T c0 is very small. This agrees with our 
instability criterion ( 1 ) that only perturbations of finite amplitude 
can grow to observable condensations (Section 5.3). 

In summary, we find that the analytic model predicts large velocities 
at the boundary between a condensation and the hot plasma in a loop. The 
growth rates of the temperature differences and the velocity are given by 
equations (5.30) and (5.45), and are plotted m Figures 5.5 and 5 . 6 . For 
fairly small perturbation amplitude, ~ the velocity becomes comparable 
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to the sound speed quite early in. the cooling process ~ before the 
condensation temperature decreases by an order of magnitude. Therefore, 
our full computer code is necessary to investigate the behavior of the 
condensation at low temperatures. 


5.4.3 Evolution of a Condensation 

As discussed m the previous section, large velocities soon build 
up due to the non-uniform cooling. With our code we can follow the 
details of the temperature, density, and velocity profiles. In this 
section we investigate the evolution of an unstable loop throughout 
its radiative cooling stage. Although we only discuss the results for 
a particular loop m this section, the general features of the evolu- 
tion turn out to be the same for all unstable loops. The only sensitive 
adjustable parameter is the amplitude of the initial perturbation; we 
describe its effect on the evolution m the following section. 

The temperature and density profiles for a typical unstable loop 
at different times of its evolution are shown m Figures 5*8 and 5*9« 

The initial profiles were assumed to be of the form 



1 + .05 cos l^-^j > (5.49b) 

where x is the Lagrangean coordinate and ± arq the positions of 

t 

the foot-points of our model. Equations ( 5 . 49a) and (5.49b) imply that 
the perturbation is approximately isobaric, centered at the top of the 
loop (x = 0), and of wavelength equal to the loop length. Since the 


n( 0 ,x) = n( 0 , 0 ) 
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conduction damping rate varies inversely as the square of the wavelength, 
equation (5.4), we expect perturbations of the maximum possible wave- 
length to be most susceptible to thermal instability. The initial values 
were- T(0,0) = 4 X 10^ K and n(0,0) = 5 X 10 10 cm" 3 . The loop was 
taken to be that due to a dipole field with a compression factor of 
T = 10 and a total length of 80,000 km, which implies a maximum height 
above the chromosphere of 30,000 km. These parameters were chosen so 
that T(0,0) ps T c , the initial temperature for the onset of instability 
(Section 5.4.1). 

From Figure 5.8 we note that by t = 480 sec the amplitude of the 
temperature variations has grown significantly, the temperature of the 
hot plasma has dropped to 1/3 of its original value, and the ratio of 
cold to hot temperatures, ? = .37- However, the ratio of minimum to 
maximum density (Figure 5-9) this time is only .73; hence, the 
temperature perturbation has begun to run away - from the density variations. 
The plasma velocities are still small, the maximum velocity is - 1.5 X 
10 cm s (directed up the loop) which implies a Mach number, M = 10" . 

We find that the condensation now enters a very rapid cooling stage, 
so that | decreases drastically. This behavior is in agreement with our 
analytical calculations in the previous section. By t = 490 sec the top 
5,000 km of the loop have dropped to chromospheric temperature, 3 X 10^ K. 
At this time the loop top becomes visible in HQ'. Comparing the temper- 
ature profile at t = 480 and t = 490 sec, we see that it requires 

R 

< 10 sec for the condensation to cool from ~ 5 X Kr K to its minimum 
temperature. (We fix the temperature so that it cannot drop below 
3 X IQ* K — Section 3*3 )• The rapid cooling is due to the strong 
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LOOP BASE 


Figure 5-9 DENSITY PROFILE OF AN "UNSTABLE’* LOOP 
AT VARIOUS TIMES IN ITS EVOLUTION. (For clarity 
the profile at t = lj-80 sec. is indicated by a 
broken line.) 





radiative losses at low temperatures; the Cox and Tucker loss function 

peaks at about 2 X IoP K. For T = 5 X 10^ K and n = 5 X 10 ® cm the 

theoretical cooling time t = ~x — ^ 1 sec. 

n A 

On the other hand, the density does not change greatly between the 
times t = 480 to t = 490 sec. Figure 5 . 9 . Hence, large pressure gradients 
result. This is to be expected. since the velocities are too small to 

move a large mass of plasma in only 10 sec. The analytic model predicts 

_2 

that for these small values of ^ 10 , the velocities should be 

highly supersonic; however, that model neglects inertial effects while 
our computer calculations indicate that they are important. The maximum 
velocity at t = 490 sec is - 4.2 X 10^ cm/sec, which is only slightly 
supersonic, but it is a factor of 3 increase over the value at t = 480 sec. 
Therefore, while the evolution of the loop is far from isobanc, it is 
also not isochoric, instead, we find that the dynamics of the plasma 
is critical m determining whether a long-lived condensation can result. 

4 

As the plasma continues to cool, the region at 3 X 10 K continues 
to expand until the complete loop reaches this temperature. This occurs 
at t = 530 sec. We note that it takes about 90 sec from the time that 

4 

the top of the loop first drops to 3 X 10 K until all the plasma cools 
to this value. If the loop were of constant density, then m HQf the 
evolution of a condensation would appear to follow the temperature 
evolution. Since the length of the loop is 80,000 km, this would 
imply an average expansion velocity of the condensation at 450 km/sec, 
which is larger than observed velocities (Section 2.2.4). However, this 
is not the case because the density profile does not remain uniform, 
instead, large variations develop as a result of the non-uniform cooling. 
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We find that the maximum upward velocity during this phase of the 

Q 

cooling occurs at t = 520 sec at a point 9 X 10° cm from, the top and is 

equal to - 7*1 X 10^ cm/sec. Prom Figure 5*8> ve can see that the 

position of this velocity is at the boundary of the condensation. The 

speed of sound' m the cold region is 3 X 10^ cm/sec; hence, a weak 

shock occurs at the boundary, M = 2. The boundary density is approximately 

10 - 3 

4 X 10 cm ; therefore, the mechanical energy flux into the condensation 
6 “2 “X 

is of order 10 ergs cm sec . However, the radiative loss rate for 

2 “1 -3 -1 

the densities m the condensation is so high, n A ^ 10 ergs cm J sec , 

that the mechanical energy flux is easily dissipated without significant 

effect on the temperature of the condensation. 

In Figure 5*9» we show the density profile at t = 570 sec, when the 

4 

temperature profile is uniform at 3 X 10 K. As a result of the dynamics, 

a large density gradient has developed at the top of the loop. The 

10 -3 

density over most of the loop volume is approximately 3 X 10 cm 

g 

however, a small, < 10 cm, clump of plasma has formed at the top with 
densities an order of magnitude higher. The clump contains 25$ of the 
total mass of loop plasma. Since the radiative losses vary as n , the 

f 

compression at the top would appear to be as much as 2 orders of magnitude 
brighter in Ha than the sides of the loop (assuming the loop is observed 
at the solar limb). 

We believe that this compression corresponds to the bright knots m 
Ha observed at the beginning of a loop prominence event. Section 2.2.4. 

The time scale for temperature variations, ^ 90 sec, is too short to 
account for the long-lived knots, ^ 10 min. But, the time scale for the 
density variations does turn out to be of the order of minutes. Section 
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5-5. Hence, we find that even though the Ha bright knots are a result 
of the thermal instability of coronal plasma, they are not a direct 
manifestation of this process. Instead, the knots are mainly due to 
the compression of loop plasma by velocities generated during the 
non-uniform cooling. As a result, their size scale is not set by the 
wavelength of the initial perturbation, which is of the order of the 
length of the loop, but is determined by the details of the plasma 
dynamics . 

The results described in this section are relatively insensitive 
to the particular parameters of a loop. The mam criterion that must 
be satisfied by a specific perturbation in order for it to develop 
into an observable knot is that its amplitude at the onset of instability 
be ^ 5 $. The importance of the plasma dynamics follows from the fact 
that the cooling time is much less than the speed of sound travel time 
across the loop at low temperatures. For the plasma m a condensation 
to be compressed significantly, the velocity generated by the cooling 
must be supersonic. If it is subsonic then, of course, the plasma 
streaming energy is less than the thermal energy; hence large pressure 
(density) enhancements will not form. 

We can estimate the dependence of the velocities on the density and 

wavelength of the initial perturbation. The velocity generated by 

2 

pressure gradients of order p/X, acting on a time scale, p/n A, is 

2 

v ~ • P . (5.50) 

m^n 0 AX 

Letting T c be the temperature at instability onset, equations (5. 50) 
and (5*9) yield 
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( 5 . 51 ) 


v « (n\) 1/7 A~ 3/7 

£ 

In the temperature range of interest, T ~ 10 K, we can approximate, 

A T \ hence, 

v « (n\) 1/3 . (5.52) 

Equation (5.52) implies that the velocity generated by the non-uniform 
cooling and, therefore, the formation of a condensation is only weakly 
dependent on n or X. Our computer calculations confirm this result. 

5 . 4.4 Effect of Perturbation Amplitude 

In agreement with the analytic model, we find that the evolution 
of the loop plasma is strongly dependent on the amplitude of the 
thermal perturbation at the onset of instability, T c « If the amplitude 
is too small, then large temperature differences do not last for a 
long enough time to produce significant density variations. 

This is illustrated by Figures 5*10 and 5*11 in which are plotted 
the temperature and density profiles of a loop with the same parameters 
as the ones used in the previous section, except that the initial 
amplitude has been halved to 2.5%. Large temperature differences occur 

4 

m this case as well, at t = 5 10 sec, the loop top is at 3 X 10 K while 

the sides are at 10^ K. However, the temperature variations last for 

less than 40 sec. By t = 550 sec, the temperature throughout the loop 
k 

is at 3 X 10 K. 

Comparing Figures 5*8 and 5*10, we note that similar temperature 
profiles occur in the loop with the 5$ amplitude and in that with the 
2,5% amplitude. But the density profiles at the time when all the 
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LOOP TOP S(l0 8 cm) LOOP BASE 

Figure 5.11 DENSITY PROFILE OF A "STABLE" LOOP 
AT VARIOUS TIMES IN ITS EVOLUTION 


ij. 

plasma m each loop reaches 3 X 10 K are very different. Comparing 
Figure 5*9 and 5«H> we note that the density profile, for the loop 
with the 2.5$ perturbation amplitude, at t = 550 sec is much smoother 
than the profile at the corresponding time, t = 570 sec, in the case 
of the 5$ amplitude. For the loop with the smaller amplitude, the 
density at the top is less than twice that at the sides and the 
compressed region is of much larger size scale. This loop would not 
appear to have a very bright, compact knot in Ha; it would appear to 
be of approximately uniform brightness instead. It would correspond 
to a loop prominence event that does not begin as a bright knot. 

For initial amplitudes larger than 5$, the evolution of the loop 
does not change qualitatively from the 5$ case. The knot appears 
brighter, smaller and lasts longer, but the overall behavior is the 
same. For example, assuming a perturbation amplitude at instability 
onset of 10$ in a loop with the same -parameters as before, results m 
relatively larger density ratios between the compressed region and 
the sides, 50 instead of 10). In this case the temperature variations 
last for 180 sec, which is a significant fraction of the lifetime of 
knots . 

The maximum velocities m the 10$ case occur at the condensation 
boundary as m the previous cases. They attain a value of -1.4X 10^ 
cm/sec, implying a Mach number of M = 4. For the case of a 5$ ampli- 
tude, we obtained M = 2, and for the 2 . 5 $ case, the maximum velocity 
is 3*2 X 10 cm/sec, i.e. Mrs 1. Hence, the velocity is approximately 
proportional to the perturbation amplitude (for small amplitudes) as 
expected from the analytic model, equation (5.48). 
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5.5 Post-Cooling Stage 

We define the third stage of the evolution of a oost-flare loop 

to be that phase during which the temperature remains approximately 

h 

constant at 3 X 10 K, Section 5.1. This evolution is, of course, 

determined by the pressure and velocity profiles that are present at 

the end of the radiative cooling stage and by gravity. We find that 

there are two possible general forms that the third stage evolution can 

take, depending on whether thermal instabilities have occurred. 

If the loop is stable, i.e. no condensations form, then as noted 

m the previous section the density profile is approximately uniform 

at the time that the cooling ends. In addition, the velocities generated 

by the cooling process are small, « c. Hence, the subsequent behavior 

of the plasma is dominated by gravity. Since the gravitational scale 
kT 

height, m~g~ s at Promos pheric temperatures is much less than the size 
of typical loops, the loop plasma begins to_fall freely. The exact motion 
is complicated due to the variation of the gravitational force along the 
length of the loop; however, as one would expect, we find that the 
maximum velocities occur at the base of the loop. They are of the order 

of the free-fall velocity, and the time scale for their growth is of 

q 

the order of the free-fall time. For a loop of height, H = 3 X 10 cm, 

the free-fall velocity is v f = 10 cm/sec, and the free- 

fall time is U/vr- « 300 sec, 
x 

If, on the other hand, the loop is unstable and a condensation 
forms, then we find that large pressure gradients are created which 
have an important effect on the velocity profile. Supersonic velocities 
form near the top of the loop well before gravitational effects become 
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significant, and at heights where the free-fall velocity is small. They 
are a result of the rapid expansion of the compressed region. In this 
section we describe m detail the post-cooling evolution of the unstable 
loop discussed in Section 5 .4.3 and compare it with the evolution of 
a stable loop. 

The radiative cooling stage ends at t = 570 sec for the unstable 
loop with the 5$ perturbation amplitude. Section 5»^-*3. The density 
profile at this time is shown in Figure 5*9 and the velocity profile 
in Figure 5*12. The velocities are residual from the non-uniform 
cooling stage; hence, within 25»000 km of the top they are directed 
up the loop. Near the base they are directed downward due to material 
having condensed onto the chromosphere. The velocity profile indicates 
that the condensation continues to be compressed (and brighten m HQ 1 ) 
for some time after the cooling ends. 

Figure 5-9 implies the existence-of- large pressure gradients which 
act to reverse the upward velocities. Since the temperature is uniform, 
the pressure profile is, of course, directly proportional to the density. 
The time required to reverse the velocity is of order 


T 


£vl 

P 


(5.53) 


where X is the width of the compressed region, £ 10 cm in our case. 
From Figure 5*12, the upward velocity is approximately - 6 X 10^ cm/sec; 
hence, M « 2 and equation (5-53) yields, T = 60 sec. This time is 
shorter than the time scale over which gravity is effective. 

In agreement with this estimate for t, we find that the velocities 
reverse direction at t = 620 sec (50 sec after the cooling ends). The 
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v (cm /sec) 



Figure 5.12 VELOCITY PROFILES OF A "STABLE" 

AND AN "UNSTABLE" LOOP AT VARIOUS TIMES IN 
THEIR EVOLUTION The broken lines refer 
to the velocity profile of the unstable loop 
at the onset of the post-cooling stage, t = 570 
sec. The broken line with no dots indicates 
negative velocities, directed uj> the loop, the 
line with dots indicates positive velocities, 
directed down the loop. 
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11 -3 

compression reaches its maximum density, 6 X 10 cm and its minimum 
7 

width, ^ 5 X 10 cm, at this time. Figure 5. 13. We note that the 

10 -3 

minimum density m the loop is approximately 2 X 10 cm at a position 

Q 

1.3 * ICr cm from the top. Therefore, the brightness m Ha of the 
compressed region should be almost 3 orders of magnitude larger than 
that of the sides of the loop; however, since the condensation is 
optically thick in H», radiative transfer effects limit its brightness. 

We also plot the density profile at t = 620 sec as a function of 
Lagrangean coordinate, x. Figure 5.14. Note that this profile is 
much smoother than the corresponding one m Figure 5*13 — the compressed 
region is clearly visible and is of almost uniform density. Over 

30$ of the loop plasma is m the condensation; however, it occupies 

-2 

less than 10 of the loop volume. Figures 5*13 and 5.1^- clearly 

illustrate the necessity of using Lagrangean coordinates for our 

computer calculations. If we were to_at tempt to follow the density 

m Eulerian coordinates, we would need a grid spacing finer than 10^ cm, 

which implies using almost 10 points over the length of the loop. In 

fact, we only use 64 points with the Lagrangean grid. 

The density profiles of the loop at t = 620, 720, 840 and 960 sec 

are plotted in Figure 5.13* The rapid expansion of the condensation 

can clearly be seen; the velocity of the condensation boundary is 

£ 

approximately 5 * 10 cm/ sec which is supersonic. 

We plot the velocity profile at t = 840 m Figure 5*12 and, for 

comparison, we also plot the profile at t = 840 sec for a loop with 

the same geometry but with no initial perturbation, i.e. a completely 

q 

stable loop. Below 10 cm from the top, the velocities in the stable 
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Figure 5.13 DENSITY PROFILE OF AN "UNSTABLE” LOOP 
AT VARIOUS TIMES IN THE POST-COOLING PHASE OF 
ITS EVOLUTION. 


90 


n \cm 



Figure 5.14 DENSITY IN AN "UNSTABLE" LOOP AS 
A FUNCTION OF LAGRANGEAN COORDINATE, x, 

AT TEE TIME WHEN THE CONDENSATION DENSITY 
IS MAXIMUM. 




loop are larger than the ones m the unstable case. This is because 
no upward velocities are generated during the radiative cooling stage 
m the stable case. On the other hand, within 10^ cm of the top, the 
velocities m the unstable loop are higher due to the expansion. The 

Q 

maximum expansion velocity at s = 8 X 10 u cm is almost equal to the 

8 

maximum free-fall velocity at s = 32 X 10° cm for the unstable loop. 

If viewed in an HQ? movie, these two loops would appear to have quite 

a different behavior. In the stable loop material would appear to 

accelerate uniformly as it fell down the sides, while in the unstable 

loop, the plasma would seem to explode out of the HQ? bright knot, 

reaching a large fraction of its maximum velocity near the top. 

Once the condensation expands sufficiently that pressure gradients 

become negligible, the subsequent evolution of the plasma is dominated 

by gravity. For the unstable loop, we find that the pressure (and 

hence the density) becomes monotomcally decreasing with height by 

10 -3 

t = 20 min. The density at the top is only-1.9 X 10 cm at this 

10 -3 

time. By t = 23 mm, the top density falls below 10 cm and we 
terminate our calculation. For the stable loop, this density is reached 
sooner at t = 18 min, however, the maximum free-fall velocity in both 

■7 

cases is approximately 10 1 cm/sec. 

In summary, we find that the main difference in the evolution of 

a stable and unstable loop is m the appearance of bright knots and in 

the acceleration of cold plasma near the top. Since we somewhat 

4 

arbitrarily fixed the temperature at 3 X 10 K, the details of our 
results on the post-cooling phase of flare loops may be inaccurate, 
but we believe that the important features are correct. We compare 
these features with observations m the following section. 
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6. DISCUSSION 


6.1 Comparison with Observations 

If we compare the results m the previous chapter on the evolution 

of an unstable loop with the description of loop prominences given m 

chapter 2, we note that they are in good agreement. For the case of the 

loop with a 5$ initial perturbation amplitude, the average density in the 

11 -3 

condensation during its lifetime is ^10 cm , and its size is small 

g 

compared to the size of the loop, ^ 5 X 10 cm. In Ha it would appear 
as a bright knot. The plasma velocities near the condensation are large, 
there is a shock at the boundary during its formation and a large 
expansion velocity during its decay. This agrees with Jefferies and 
Orrall's ( 1965 a) observation of large streaming velocities at the top 
of active loop prominences. Defining the lifetime of a knot as the 
time interval from the instant it first appears in Ha to when it expands 
to a size of half the loop length; we~find that the lifetime in our 
case is approximately 8 minutes. Bruzek and Kuperus (1972) quote the 
lifetime of a knot as 3 “ 10 minutes. 

The maximum velocity that we obtain at the loop base, 115 km/sec, 
equals the free-fall velocity. The average velocity is slightly lower, 

^ 100 km/sec, in agreement with observations. There is a shock at the 
base where the falling matter impacts with the chromosphere; however, 

the radiative losses are sufficient to dissipate the mechanical energy 

* 

flux. In an unstable loop, large velocities appear near the top (at 
the knot boundary) due to pressure gradients, whereas, in a stable loop 
material accelerates uniformly down the sides. We calculate the total 
lifetime of a loop event as approximately 30 minutes. Bruzek (1964a) 
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gives the observed lifetime of typical loops as 30 - 60 minutes . Of 
course, these time scales and velocities will vary depending on the 
size and density of a loop, equations (5*2) and (5-3 )• For lower 
densities and larger loop heights, the radiative cooling time and 
the free-fall time increase. However, the long lifetime of loop 

5 

prominence systems, <S 10 sec, does not appear to be due to a cooling 

effect. Even if all conductive losses are suppressed by a large 

magnetic field compression factor T, the radiative losses alone for 
10 -3 7 

densities of ~ 10 cm and temperatures of 10 K imply a cooling 

3 

time of only 10 sec. Hence, the long lifetime of a loop system is 
most likely due to a long-lived flare heating process. 

In general we find that the temperature, densities and velocities 
measured m a loop prominence are accurately reproduced by our model. 

From this, we conclude that the thermal instability mechanism by itself 
is sufficient to account for a loop prominence event. The appearance of 
knots and their subsequent evolution can be understood as being due to 
' the temperature and density dependence of the cooling function (Figure l.l) 
for the coronal plasma. Our only important assumption is the form of 
the temperature and density profiles during the radiative cooling stage. 

If a perturbation of small but finite amplitude is present a condensation 
will form. If, on the other hand, the profiles are uniform as we would 
expect from the discussion in Section 5*2.1, then no knots appear and 
the evolution of the loop plasma will more closely resemble that of 
coronal rain. 

We have given no justification for choosing the particular form 
of the perturbation m equation (5»^9)» except that it is the one with 
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the largest wavelength and, hence, the fastest growth rate. In accordance 
with observations, condensations may form at the sides of loops if the 
perturbation profiles are suitably chosen; however, these knots are not 
as dense or as long-lived as those centered at the top. Unfortunately, 
the observations are not accurate enough to determine the temperature 
and density profiles of a loop to within hence we cannot directly 
check our assumption. The soft x-ray emission is observed to be 
greater at the tops of loops (Cheng and Widmg, 1975, Pallavacini et al., 
1975 )> hut that would be expected due to the area factor even if the 
temperature were slightly lower there. 

There is reason to expect that the temperature will be higher at 
the sides of a flare loop if the heating is due to high energy electrons 
as m the models proposed by Sturrock ( 1968 ), Petrosian (1973) and others 
(Svestka, 1976, Section 3*2). In these models the bulk of the flare 
energy is carried down the field lines by non-thermal electrons and 
deposited in the chromosphere or lower corona. Hence, the temperature 
should initially be higher near the base of a loop. If the electrons 
have non-zero pitch angles and are trapped m the flux tube, then 
again, we would expect higher temperatures at the sides of a loop. The 
electrons spend more time near their mirror points than near the top of 
a loop, since their velocity parallel to the field is highest at the top. 
This process also has the attraction that if a loop is symmetrical about 
the top, then the heating due to trapped energetic electrons will have 
the same symmetry; hence any temperature perturbation will be symmet- 
rical about the top m agreement with our assumption m equation ( 5 *^- 9 ) • 

We intend to investigate these possibilities in detail by including 
a heating function in our code to represent flare heating. We also 
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intend to modify the boundary conditions to permit an upward mass flux 
at the base of the loop. This will allow us to study the evaporation 
of chromospheric material, and, perhaps, to determine the form of 
flare heating. 

6.2 Conclusions 

To summarize, we restate the main conclusions of our calculations 

here 

1. For typical post-flare loops, the conductive damping time at high 

7 

temperature, ^10 K, is so short that all temperature and pressure 
gradients are damped quickly compared to the cooling time. Hence the 
analytic models proposed by Antiochos and Sturrock (1976a, b) adequately 
represent the initial cooling. 

2. The inclusion, of the area factor does not alter the linear growth 
rate of thermal perturbations, equation (l.l4). 

3. The effect of the nonlmearities is to_increase the stability of 
the plasma. Figure 5»3- 

4. Observable condensations will grow if criteria 1 - 4 m Section 5»3 
are satisfied. The important requirement is that the perturbation, 
amplitude at the onset of instability be ^ 5$* 

5. Shocks form at the boundary of knots during their compression phase. 
Large velocities are generated by their subsequent expansion. 

6. The appearance of bright knots is mostly due to density variations 
along the loop, since temperature variations last for only a short time, 
< 100 sec. 

7- The evolution of loop prominences can be understood as simply due 
to thermal instability in chromospheric material that has been 
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evaporated into the corona by a flare. Mechanisms such as strong 
compression by magnetic fields or the migration of energetic protons 
across field lines are not required. 

8. The radiative cooling time for temperatures and densities typical of 
post-flare loops is too short to account for the long lifetime of loop' 
prominence systems; hence, a long-lived (time scales of the order of 
several hours), flare heating process is required. " 
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APPENDIX A 


FLUID EQUATIONS FOR A DIPOLE LOOP 


We wish to derive the fluid equations for the case of a flux tube. 
Since we assume that the field dominates the plasma, we can most easily 
incorporate the effects of the field by letting the field lines define 
a curvilinear coordinate system. This is always possible if the field 
is potential (Morse and Feschback, 1953, P. 14). 

Let (x , x , x J ) denote the coordinate system, where x = s, the 
distance along the field. We are only interested m an individual 
loop whose cross-sectional area is very small compared to the radius 
of curvature of the field. Thus, we can approximate the metric in 
the loop by 



(Al) 


To obtain the fluid equations in terms of our coordinate system, we first 
rewrite them in covariant form. They are* 


(1 i 

all “bV 

St 


■ ( A2) 

If™ 11 1) + Kv 1 + p 6 i + n ^ p ■ “"h®! * (A3) 

+ I p) + (| + | pv^ + q J + nv) 

' llJ = run^v + £ 


(A4) 


where we have used the form of the equations given by Braginskn ( 1965 ) 
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and the tensor notation of Adler, Bazin and Schiffer ( 1965 ), Chapter 2. 
Note that all spatial derivatives are covariant derivatives and that the 
Einstein convention of summing over repeated indices is used. 

The equations A2 - A4 simplify considerably in our case since we have 

2 3 

no variables dependent on x or x and all vectors have a component in 

1 • 
the x direction only. Hence, 

v 1 = (v (x 1 ), 0, 0 ) . ( A 5) 


The heat flux and the stress tensor H are given by Braginskn 

(1965); 


and 





(A6) 


n 

ij 


n 




j-iij 


+ V 


Jill 


2 

3 8 ij 



(A7) 



Using equation (A 5 ) and the metric (Al), the mass equation (A2) becomes 


3n 1 aiAnv) _ 

St A ds 


(A9) 


where we have replaced x^" by s. Using (A 9 ) and (A8), the momentum 
equation (A3) becomes 
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+ 


v3v \ S£ 1 3_(Ar|w ) 

3s/ 3s A 3 s 


1 3 (r)w ) 

" 3 37 = V®, 


(A10) 


where 


23v v 3 An A 
W 3s 3s 


(All) 


and we have replaced by g . 

Finally, using equations (A10), (A 9 ) and (A 6 ), the heat equation 
(A4) can be expressed as 


^/3£ v3p.\ ^£ 3_(Av) _ 1 3_(K3t) X] w 2 

2 \3t 3s / 2 A 9s A 3s 3s 3 


(A12) 


Since T| is small, we can neglect the viscosity terms to obtain the set 
of equations (3.20) - (3-22) in Section 3-6. 

We now need to calculate A(s) and g^(s) f° r a dipole field. Using 
the coordinates defined in Figure (3.1), a dipole field is given by 
(Jackson, 1962 , p. 143), 

B ec (2 sm 0 - cos 0 1 Q ) . (A13) 

The equation of a field line of height R is given by 

r/R = cos 2 0 . (Al4) 


Letting 


NJ = sm 0 , 


(A15) 


the distance along the field line, 

ds =T^( dr ) 2 + r 2 (de ) 2 , 


(A16) 
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becomes * 


ds 


= *\f: 


3v *t* 1 dv . 


(A17) 


Integrating (A17), we obtain equation (3.28). The area factor equation 
(3.29) can be obtained from equations (A13) and (Al4) since we have 


A cc 


PI 


(A18) 


The component of gravity parallel to the field, equation (3 .30), is 
given by 


“ ' g O 


l . 


B 

PI 


(A19) 


where JL is a unit vector in the vertical direction, 
2 


i = cos 0 H - sin 6 I 
z r 0 


(A20) 


In Figure Al, we plot A(s) and g (s) for a dipole loop. 
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Figure A. 1 PROFILES OF THE AREA OF A LOOP DUE TO A 
DIPOLE SOURCE AND OF THE COMPONENT OF GRAVITY PARALLEL 
TO THE LOOP. The solid line refers to A(S) and the 
broken line to g. (S). Both A and g„ are normalized to 
equal unity at their respective maxima. 
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APPENDIX B 


C B.l FINITE DIFFERENCE EQUATIONS 
C 

C DETAILS OF DIFFERENCING SCHEME FOR EQUATIONS . 3 - 4 . 7 ) WITH 
C BOUNDARY CONDITIONS ( 4.18 - 4 . 19 ). THE FOLLOWING IS THE CODE TO 
C ADVANCE VALUES OF VARIABLES BY ONE TIME STEP, DT, THE VALUE OF DT 
C FOR THE NEXT TIME STEP CAN BE EVALUATED USING CONDITION ( 4 . 40 ). 

C THAT ROUTINE IS NOT INCLUDED HERE. IN THE NEXT SECTION A ROUTINE 
C FOR TABLE LOOK-UPS IS GIVEN. 

C 

C WE DIVIDE EACH LEG OF THE LOOP INTO 32 EQUAL MASS SLABS. EACH SLAB 
C IS OF WIDTH (MASS) DX. WE CHOOSE THE POINT X = 0 TO BE THE TOP OF 
C THE LOOP. 

C 

C LIST OF VARIABLES USED IN CODE: 

C ALL ARRAYS HAVE 33 ELEMENTS. (SOMETIMES NOT ALL ELEMENTS USED) 

C PRIMARY VARIABLES' 

C VALUES OF PLASMA PARAMETERS AT TIME = T. NEEDED AS INPUT TO ADVANCE 
C TIME STEP. MUST BE OUTPUTTED TO USE FOR THE NEXT STEP. 

C "S(J)" = EULERIAN POSITION OF LAGRANGEAN POINT (j-l)*DX 
C "V(J) n = VELOCITY AT LAGRANGEAN POINT (j-l)*DX. 

C "T(J)" = TEMPERATURE AT LAGRANGEAN POINT (J-1/2)*DX. 

C (NOTE THAT WE POSITION SOME VARIABLES DIFFERENTLY FOR PROPER 
C CENTERING OF DIFFERENCE EQUATIONS.) 

C "P(J)" = PRESSURE AT LAGRANGEAN POINT (j-l/2)*DX. 

C "AN(J )" = NUMBER DENSITY AT LAGRANGEAN POINT (J-1/2)*DX. 

C "Q (j )” = ARTIFICIAL VISCOSITY AT LAGRANGEAN POINT (j-r/2)*DX. 

C 

C PROVISIONAL VARIABLES 

C VALUES OF PLASMA PARAMETERS AT TIME = T + DT/2. THE SPATIAL 
C POSITIONING IS SAME AS ABLVE. 

C "S1(J)", "V1(J)", "T1(J)'\ "Pl(j)", "AN1 (J )" , "Q1 (J )" 

C 

C SECONDARY VARIABLES 

C VALUES OF SECONDARY PARAMETERS AT EITHER TIME = T OR T + DT/2. 

C MAY BE OUTPUTTED IF DESIRED. 

C "VDER(j )" = D(V(J))/DT. 

C "TDER(j)" = D(T(J))/DT. 

C "AS (j )" = AREA FACTOR AT LAGRANGEAN POINT (j-l)*DX 
C "GS (J )" = COMPONENT OF GRAVITY PARALLEL TO 
C LOOP AT LAGRANGEAN POINT (j-l)*DX. 

C "T12(J)" = SQUARE ROOT OF T(j). 

C "T72 (J )" = T(J) RAISED TO THE POWER OF 7/2. 

C "ALT(j )" = COX AND TUCKER RADIATIVE LOSSES AT TEMPERATURE T(j). 

C "HFLX(J)" = HEAT FLUX AT LAGRANGEAN POINT (J-l)*DX. 

C 

C ARRAYS USED FOR STORAGE ONLY "V12(j)'\ ”AV(j)", "D(j)". 

C 

C CONSTANTS CHOSEN TO MAKE VARIABLES DIMENSIONLESS. 

C "DX" = 2/3, "EQV1" , "QOPO", "EQT3", "EQT32". 

C "AIB" = VOLUME OF LEG OF LOOP 

C "DT" = TIME INCREMENT, CONSTANT ONLY FOR EACH TIME ADVANCE. 

C "DT2" = DT/2, "DT3" = DT/3, "DT4" = DT/4. 
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C VARIABLES USED FOR STORAGE ONLY' "AIL, AIF , DELTV" 

C 

C BOUNDARY CONDITIONS VALUES CONSTANT FOR ALL TIME 

c s ^ j? ’ J ^ j x ^LX(1 ) , (33 ) , 3 V(33 ) , J(3 3 ) _ ^ 

C ROUTINE TO INCREMENT TIME BY DT/2. 

C CALCULATE PROVISIONAL VARIABLES* VI (2), SI (2), ANl(l), Ql(l) 

VDER(2) = EQV1*AS(2)*(P(2) + Q(2) - P(l) - Q(l)) + GS(2) 

VI (2-) = V(2") + DT2 x_ VDER(2 ) 

V12(2) = Vl(2) + V(2 ) 

SI (2) = S(2) + DT4*V12(2) 

C— INSERT ROUTINE TO CALCULATE AS (2), GS(2), AIL - SEE SECTION B.2 
C AIL = VOLUME OF LOOP FROM TOP TO POINT SI (2) 

310 AN1(1) = DX/(AIL) 

AV(2) = AS (2)*V1 (2) 

DELTV = VI (2) - VI (1) 

IF (DELTV )330, 320, 320 
320 Ql(l) = 0. 

GO TO 340 

330 Q1 (1 ) = QOPO*ANl(l)*DELTV*DELTV 
340 CONTINUE 

C START MAIN DO LOOP 111111111111111111111111111111111111111 
DO 460 J = 1, 30 

C CALCULATE Vl(J+2), Sl(J+2), ANl(j+l), Ql(j+l) 

VDER(j+2) = EQVl*AS(j+2)*(P(J+2) + Q(J+2) - P(j+l) - Q(J+l)) 

1 + GS ( J+2 ) 

Vl(j+2) = V (J+2 ) + DT2*VDER(J+2) 

V12 ( J+2 ) = VI (J+2) + V ( J+2 ) 

SI (J+2) = S ( J+2 ) + DTA*V12 ( J+2 ) 

C INSERT ROUTINE TO CALCULATE AS(j+2), GS(J+2), AIF - SEE SECTION B.2 

C AIF = VOLUME OF LOOP FROM TOP TO POINT SI ( J+2) 

370 ANl (J+l) = DX/ (AIF- AIL) 

AIL = AIF 

AV(J+2) = AS (j+2 )*V1 (J+2 ) 

DELTV = VI (J+2) - VI (J+l) 

IF (DELTV) 390, 380, 380 
380 Q1 (J+l ) = 0. 

GO TO 400 

390 A1 ( J+l ) = QOPO^ANl ( J+l )*DELTV*DELTV 
400 CONTINUE 
C CALCULATE Tl(j), Pl(J) 

D (j+l ) = AS(j+l)*AS(J+l)*ANl(j)*ANl(j+l)/(ANl(j) + ANl(j+l)) 

KFLX (J+l ) = D(J+1)*(T72(J+1) - T72(J)) 

TDER(J) = (P(J) + Ql(j))*(AV(j) - AV ( J+l ) ) + ANl (j)*ALT(j) 

1 + EQT3*(HFLX(j+l) - HFLX(j)) 

T1(J) = T(j) + DT2*TDER(J ) 

IF (T1 (J ) .LT. T(33))GO TO 440 

C INSERT ROUTINE TO CALCULATE T12(j), T72(j), ALT(j) - SEE SECTION B 2 

GO TO 450 

440 Tl(j) = T(33) 
alt(j) = alt(33) 

T72(J) = T72 (33 ) 

450 P1(J) = ANl (J )* (T1 (j) ) 

460 CONTINUE 

C END OF MAIN DO LOOP 111111111111111111111111111111111111111 
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C CALCULATE AN1(32), Ql(32) 

AN1 (32) = DX/ (AIB-AIL) 

DELTV = VI (33) - VI (32) 

IF(DELTV) 480, 470, 470 
470 Q1 (32) = 0. 

GO TO 490 

480 Ql(32) = Q0P0*AN1(32)*DELTV*DELTV 

490 CONTINUE 
C CALCULATE Tl(3l), Pl(3l) 

D(32) = AS(32)*AS(32)*AN1(31)*ANI(32)/(AN1(32) +ANl(31)) 

HFLX(32) = D(32)*(T72(32) - T72(31)) 

TDER(3I) = (P(31) + Q1(31))*(AV(31) - AV(32)) + AN1 (31 )"*ALT (31 ) 

1 + EQT3* (HFLX(32) - HFLX(3l)) 

TI(31) = T(31) + DT2*TDER(3X) 

IE (T1 (31 ) .LT. T(33))G0 TO 530 

C INSERT ROUTINE TO CALCULATE T12(3l), T72(3l), ALT(3l) - SEE SECTION B 2 

GO TO 540 

530 T1(3L) = T(33) 

ALT (31) = ALT (33) 

T72(3l) = T72 (33 ) 

540 PI (31) = AN1 (31 )* (T1 (31) ) 

C CALCULATE Tl(32), PI (32) 

hflx(33) = as (33 )* (T72 (33 ) - T72(32))/(s(33) - si (32)) 

TDER(32) = (P(32) + Q1(32))*AV(32) + ANI (32 )*ALT(32) 

I + EQT32*HFLX(33) - EQT3*HFLX(32) 

Tl(32) = T(32) + DT2*TDER(32) 

IF (T1 (32) .LT. T(33)) GO TO 580 

C-— INSERT ROUTINE TO CALCULATE T12(32), T72(32), ALT (32) - SEE SECTION B 2 
GO TO 59O 

580 T1 (32 ) = T(33) 

/ ALT(32) = ALT (33 ) J 

T72(32) = T72(33) 

590 PI (32) = ANI (32 )* (T1 (32 ) ) 

C END OF TIME INCREMENT BY DT/2. BEGIN INCREMENT BY DT. 

C CALCULATE PROVISIONAL VARIABLES* V(2), S(2), AN(l), Q(l) 

VDER(2) - EQV1*AS(2)*(P1(2) + Ql(2) - Pl(l) - Ql(l)) + GS(2) 

V(2) = V(2) + DT^VDER(2 ) 

S (2) = S(2) + (V12(2) + V(2))*DT3 

C— INSERT ROUTINE TO CALCULATE AS(2), GS(2),AIL - SEE SECTION B.2 
C AIL = VOLUME OF LOOP FROM TOP TO POINT S(2) 

620 AN(l) = DX/ (AIL) 

DELTV = V(2) - V(l) 

if(deltv)64o, 630,630 
630 Q(l) = 0. 

GO TO 650 

640 Q(l) = Q0P0*AN ( 1 )*DELTV*DELTV 
650 CONTINUE 

C START OF MAIN DO LOOP 2222222222222222222222222222222222222 
DO 770 J=l, 30 

VDER(j+2) = EQVl^AS (j+2 )* (PI (j+2 ) + Ql(j+2) - Pl(J+l) - Ql(J+l)) 

1 + GS ( J+2 ) 

V (j+2 ) = V (j+2 ) + DT*VDER ( J+2 ) 

S(J+2) = S(j+2) + (V12(J+2) + V(J+2))*DT3 
C— INSERT ROUTINE TO CALCULATE AS(j+2), GS(J+2>, AIF - SEE SECTION B.2 
C AIF = VOLUME OF LOOP FROM TOP TO POINT S(j+2) 
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680 AN(J+1) = DX/(AIF-AIL) 

AIL = AIF 


DELTV = V (J+2 ) - V('J+l) 

if(deltv) 700, 690, 690 
690 Q(J+1) = 0. 

GO TO 710 

700 Q ( J+l ) = qopo*an (j+i )*deltv*deltv 

710 CONTINUE 
C CALCULATE T(j), P(j), 

.HFLX.(Jrfel-)- = D(-J+l-)*-(T72(j+l) - T72 ('J')') 

TDER(J) = (P1(J) + Q1(J))*(AV(J) - AV(J+l) ) + AN1(J)*ALT(J) 

1 + EQT3*(HFLX(j+l) - HFLX(J)) 

T (j) = T(j) + DT*TDER(j) 

IF (T.(J ) .LT. 'T(33'))GO TO 750 

C— INSERT ROUTINE TO CALCULATE T12(j), T72-(j), ALT(J') - SEE SECTION B.2 
• GO TO 760 
750 T(J) = T(33) 

T12(J) = T12(33) 

T72 ( J ) = T72(33) 

760 P(J) = AN (J )*(T(j ) ) 

770 CONTINUE 

C END OF MAIN DO LOOP 2222222222222222222222222222222222222222 
C CALCULATE AN (32), Q(32) 

AN (32 ) = DX/ (AIB-AIL) 

C CALCULATE T(3l), P(3l) 

DELTV = V(33) “ V(32) 

IF (DELTV ) 790 , 780 , 780 
780 Q(32) = 0 
GO TO 800 

790 Q (32) = QOPO*AN (32 )*DELTV*DELTV 
800 CONTINUE 


840 


850 


ALT (31 ) 
T12 (31 ) 
T72(31) 
p(31) = 


HFLX (32 ) = D(32)*(T72(32) - T72(3l)) 

TDER(31) - (PI (31) + Q1(31))*(AV(3.1) - AV(32)„) + ANl(3l)*ALT-(3l) 

I + EQT3*(HFLX(32) - HFLX (31)) 
t(31) = t(31) + dt*tder(31) 

IF (T(‘3l) .LT. T(33))G0 TO 84o 

C— -INSERT ROUTINE TO CALCULATE 112(31), T72(3l), ALT (31) - SEE SECTION B.2 
GO TO 850 
T(-31) = T(33) 

» alt ( 33) 

= T12(33) 

= T72(33) 

AN (31 )* (T (31 ) ) 

C CALCULATE T(32), P(32) 

855 HFLX (33 ) = AS (33 )*(T72-(33 ) - T72(32))/(s(33) - S1(32)J 
TDER(32) = (PI (32) + Q1 (32) )*AV (32) + ANI(32),*ALT(32), 

1 + EQT32*HFLX(33) - EQT3*HFLX(32) 

, T(32) = T(32) + DT*TDER(32) 

IF(T(32) .LT. T(33))GO TO, 89O 

C— INSERT ROUTINE TO CALCULATE T12(32), T72(32), ALT (32) - SEE SECTION B.2 
GO TO 900 
890 ’ T(32) = T(33) 

ALT (32) = ALT (33) 

T12(32) = T12(33) 

T72 (32) = T72 (33 ) - 

900 P(32) = an( 32 )*(t( 32)) 
c END OF ROUTINE * 

C ALL PRIMARY VARIABLES HAVE BEEN INCREMENTED BY DT. 

C ^x^x^^^x*^*-****^** ** X X X x x x x x x x x *x- x x - x - x ttx-x^- x -xx- x -xx 
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C B.2 ROUTINE FOR TABLE LOOK-UPS 
C 

C AS DISCUSSED IN SECTION 4.4 WE MUST EVALUATE 3 FUNCTIONS OF EULERIAN ' 

C POSITION S (j) ; AS(j), GS(j), AND AIF (OR AIL), AND 2 OF TEMPERATURE 
C T(j)* T12 (j ) AND ALT (J ) . WE FIRST TABULATE THESE 5 FUNCTIONS AT 240 
C LINEARLY SPACED VALUES IN EACH OF THE FOUR INTERVALS- 
C S(J) OR T (j) = 1/4096 - 1/256, 1/256 - 1/16,1/16 - 1, AND 1 - 16, 

C AND THEN USE LINEAR INTERPOLATION TO FIND THE VALUE OF A FUNCTION FOR 
C A PARTICULAR S(j) OR T(j). BY SUITABLY NORMALIZING S(j) AND T(j) WE 
C CAN USE THE FIRST 2 BYTES OF S(j) OR T(j) TO LOCATE THE NEAREST ENTRY IN 
C THE TABLE. FOR EXAMPLE, WE LIST BELOW THE ROUTINE FOR FINDING ALT(j) 

C AND T12 (J ) . IN THIS ROUTINE THE TABLES ARE IN ARRAYS AL AND TSQ , THE 
C POSITION OF THE NEAREST ENTRY IS EITHER n NST" OR "NST+l". 

C 

REAL*4 TCH, STREF, TABSEP 
LOGICAL*! QQ,QL(4) 
integer*4 NFSET 

INTEGER*2 NN,MM 

EQUIVALENCE (TCH, NN, QQ), (TABSEP, QL(l)), (MM, STREF) 

DATA NFSET/Z00003E0F/, TABSEP/ ZOOO lOOOO/, 

1STREF / ZOOOOOOOO/ 

C- — ROUTINE TO CALCULATE T12(j), T72(j), ALT(j) 

C FIRST WE FIND "NST", THE POSITION OF THE ENTRY NEAREST 
C TO T(J) AND LESS THAN T(j). 

TCH = T(J) 

NST = NN - NFSET 

C NOW WE FIND THE SEPARATION BETWEEN T(j) AND THE POSITION "NST". 

QL(1) = QQ 
MM = NN 

DIFST = TCH- STREF 

C NOW DECIDE WHETHER "NST" OR "NST+l" IS CLOSER TO T(j). 

C FOR EXTRA ACCURACY WE USE SECOND ORDER INTERPOLATION FOR T12(j). 

IF (DIFST+DIFST-TABSEP )720 , 720 , 730 
720 ALT(J) = AL(NST) + DIFST* (AL (NST+l) - AL (NST )) /TABSEP 
TRATO = DIFST/ (STREF+STREF ) 

T12(J) = TSQ(NST)*(1.+TRAT0*(1.-.5*TRAT0)) 

GO TO 740 

730 DIFST = DIFST-TABSEP 

ALT ( J ) = AL(NST+1) + DIFST* (AL (NST+l) - AL (NST )) /TABSEP 
TRATO = DIFST/ (STREF+STREF+TABSEP+TABSEP) 

T12 (J ) = TSQ (NST+l )*(1.+TRAT0*(1.-.5^TRAT0)) 

740 T72 ( J ) = TCH*TCH*TCH*(T12(J) 

C--END OF ROUTINE 
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